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ABSTRACT 

We present a simulation analysis of flexion and shear measurement using shapelet decompo- 
sition, and identify apparent qualitative differences between flexion and shear measurement 
noise in deep survey data. Taking models of galaxies from the Hubble Space Telescope Ultra 
Deep Field (HUDF) as a basis set and applying a correction for the HUDF PSF we gener- 
ate lensed simulations of deep, optical imaging data from Hubble 's Advanced Camera for 
Surveys (ACS) with realistic galaxy morphologies. We find that flexion and shear estimates 
differ in our measurement pipeline: whereas intrinsic galaxy shape is the largest contribution 
to noise in shear estimates, noise in flexion estimates is dominated by pixel noise due to fi- 
nite photon counts and detector read noise. This pixel noise also increases more rapidly as 
galaxy signal-to-noise decreases than is found for shear estimates. We provide simple power 
law fitting functions for this behaviour, for both flexion and shear, allowing the effect to be 
properly accounted for in future forecasts for flexion measurement. Using the simulations we 
also quantify the systematic biases of our shapelet flexion and shear measurement pipeline for 
deep Hubble datasets such as GEMS, STAGES or COSMOS. Flexion measurement biases are 
found to be significant, but consistent with previous studies. 
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1 INTRODUCTION 

The study of weak gravitational lensing has progressed to be- 
come one of the most important t echniques in observatio nal cos- 
mology (see e.g. lSchneide"rll2006l : iHoekstra & Jainl|2008l for a re- 
view). As it does not depend upon the microscopic composi- 
tion of the mass by which it is caused, gravitational lensing re- 
quires no assumptions to be made regarding baryon-dark mat- 
ter physics. It can thus be used to make dire ct observations 
of the matter distribution on large scales (e.g . IHoekstra et al.l 
20061: iMassev et all 120073: iBeniamin et alj|2007l : iFu et alj|2008l ; 
Schrabback et al. 20 id : Heymans et al. 20121 : Kilbinger et al., in 



prep; Benjamin et al., in prep.). 

Weak lensing studies have typically measured the small but 
coherent distortions in the ellipticities of distant galaxies, due to 
the shear 7, and have used these measurements to constrain the dis- 
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tribution of the intervening matter field. The weak lensing descrip- 
tion has been extended to higher order via the fl exion formalism 
(see iGoldberg & Baconll2005l : iBacon et alj|200fj . hereafter B06), 
which describes the slight arcing of galaxy shapes. Despite being 
a weaker effect than shear, it has been hoped that the flexion sig- 
nal may yet be profitably measured - galaxies typically display less 
intrinsic curvature than intrinsic ellipticity, and so the contribution 
to noise from intrinsic shape is reduc ed for measurements of flex- 
ion. S t udies have suggested ( e. g. B06;IOkura. Umetsu & Futamasel 
20071: lLeonardetal.1 120071; lOkura. Umetsu & Futamasd 120081: 



Leonard. King & Goldberg! 1201 1 ) that flexion measurements can 
provide extra constraints upon galaxy-galaxy lensing results and 
cluster mass reconstructions. 

For all weak lensing analyses, the correct treatment of system- 
atic errors is vital: image distortions and shape bias due to convolu- 
tion with an anisotropic Point Spread Function (PSF), as possessed 
by all optical instruments, are commonly an order of magnitude 
greater than the gravitational signal we wish to measure. The opti- 
mal, unbiased estimation of weak lensing signals from real data is 
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the subject of much ongoing research, involving a variety of differ- 
ent approaches towards the accurate inference about galaxy shapes, 
accounting for telescope optics, detector effects and noise. Many 
of the current methods used to corr ect for the effects of the PSF 
are ba s ed on the scheme proposed bv lKaiser. Squires & Broadhurstl 
l ll995h . lLuppino & Kaiser] dl997l) and lHoekstra et alj dl998h . com- 
monly referred to as KSB or KSB+. The use of these techniques 
has proved to be both successful and widespread, to date. 

Despite its practical success, there are certain eleme nts of 
the KSB treatment that are unsatisfactory: I Raised J2000h pro- 
vides one compelling discussion of these potential limitations. 
This has prom pted ef f orts to develop alternative weak l ensing 
methods (e.g. iKaisej |2000|: IRhodes Refregier & GrotrJ l200(f 
Bernstein & Jarvisl |2002t iRefregied 120031; IRefregier & Bacorj 
2003tlHirata et al.ll2004ljMassev & Refregierll2005l:lKuiikerj|2o6^ 
Melchior. Meneghetti & Bartelmann 2007; Nakajima & Bernstein! 



20071: iMiller et al.l |2007|; iKitching et al.l 120081; iBernsteirj 



Viola. Melchior & BartelmanrJ 1201 ll: iMelchior et al] 1201 ll : 



Kacprzak et all 120121 ; 



Miller et al 



20 id: 



20121 : Zuntz et al., in 



The Shear TEsting Program (STEP: see iHevmans et al J 1 20061: 
Massey et al] l2007bl) and GREAT challenge series dBridle et al .1 



2009 1 l20ld : IKitching et al]|201lL 120121) have compared a wide 



range of current weak shear estimation methods, using blind-tests 
on simulated lensing data. 

The s hapelet approach , proposed by iBernstein & Jarvisl 

d2002h and IRefregier! J2003I) . is one such alternative to KSB 
methods. Shapelets expresses galaxy images as a sum of sim- 
ple basis functions — Gauss-Laguerre or Gauss-Hermite poly- 
nomials — that behave well under deconvolution with a mod- 
elled PSF. In addition, the first method for the practical estima- 
tion of the flexion signal w as built within the shapel ets framework 
Goldberg & Baconl Booll : B06). Further work by iMassev et all 



2007), hereafter referred to as M07, investigated shear and 



flexion measurement with in the polar shapelets formalism of 
IMassev & Refregier (2005); results suggested that polar shapelets 
provided an apparently natural framework for e stimating both 
quantities. IVelander. Kuiiken & Schrabbackl d201ll) used shape lets 
(in an implementation closely related to that of Kuiiken 200of) to 
constrain flexion in Hubble Space Telescope (HST) data, although 
following a somewhat different modelling strategy to that sug- 
gested in M07. An alternative method for measuring flexion, re- 
ferred to as Higher Order Lensing Image Characteristics (HOLICS) 
has also been suggested, and indeed employed in cluster modelling 
from real data tofaira et al.l200ll200^ : lGoldberg & Leonard! 20071 ; 
iLeonard et al .120071 120111) to provide less noisy measurements than 
shapelets. HOLICS is conceptually a generalization of KSB meth- 
ods to higher order image moments. 

In this paper we present an analysis of simulations of 
space-based lensing data, such as that taken using the HST Ad- 
vanced Camera fo r Surveys (ACS: see, e.g.. lHartig et al.l 120031 ; 
IRhodes et af]|2007l) . Several wide-area imaging surveys that may 
be used for weak lensing exist for this instrument, including the 
Gal axy Evolution from Morphology & SE Ds survey (GEMS: see, 
e.g. jRix et al. I l2004l ; IHevmans et alJPiool), the Cosmic Evolution 
Surve y (COSMOS: see, e.g.. lScoville et al]|2007l ; lLeauthaud et al] 
120071) and the Spac e Telescope A901 /902 Galaxy Evolution Survey 
(STAGES: see, e.g jGrav et alj2009]) . The imaging data from these 
surveys share important characteristics as regards lensing measure- 
ment, having a small but non-Gaussian PSF and significant corre- 
lation in the noise between pixel s due to the dithering and driz- 
zling strategies employed (see e.g jLeauthaud et alj2007h . In addi- 
tion, there exist a wealth of galaxy cluster imaging data in the ACS 



archive that are also of interest for shear and flexion lensing anal- 
ysis. Mindful of these factors, we construct realistic simulations of 
ACS lensing data using real sky noise taken from blank areas in the 
GEMS survey data, and with a PSF that matches the radial profile 
of the GEMS PSF. Within these simulations we apply known input 
shears and flexions, and use the resulting measurements to calcu- 
late the necessary calibration for the shapelet measurement of shear 
and flexion from space. 

Our paper is organised as follows. Sections[2]and[3]begin with 
a brief description of the flexion formalism and a summary of our 
adopted shapelet measurement method. In Section [4] we describe 
how we make shapelet models of Hubble Ultra Deep Field data, 
including both stars and galaxies. These are used to generate sim- 
ulations of lensed ACS data, which we describe in Section [5] In 
Section [6] we test how well we can measure flexion on this simu- 
lated data. Finally, we present our conclusions in Section|7] 



2 WEAK SHEAR AND FLEXION FORMALISM 

To begin, we review the flexion formalism developed by 
iGoldberg & Baconl J2005T) and B06, examining briefly how weak 
flexion is defined in relation to weak shear. We restrict the discus- 
sion to lensing measurements in the we ak regime, so we do not 
consider the reduced shear or flexion (see lSchneider & Seitzl 19951 : 
ISchneider & ~Erll2008l) . 

Gravitational lensing conserves surface brightness, so the 
effects of lensing may be described in terms of coordinate 
transformations between the lensed and unlensed sky coordi- 
nate plane. In general, the relationship between these coordi- 
nates is non-linear and is de s cribed by the le ns equation (e.g., 
iBartelmann & Schneiderll200ll : ISchneiderll2006l) . 

However, if we may assume that changes in the lens proper- 
ties of a system occur only on angular scales that are large com- 
pared to the angular size of the lensed light source, then the lens 
equation may be locally linearized. In what follows, lensed and un- 
lensed coordinates are denoted by 6 and 6' respectively, and we 
define position relative to a galaxy centroid as AO = — C and 
AO' = 0' — 0' c , where C and 0' c are the coordinate centroids of 
the galaxy in the image and source plane, respectively. Approxi- 
mately linearizing the lens equation around this centroid, we write 



AO- ~ AijAOj, 



(1) 



where Aij is the Jacobian matrix of the transformation given by the 
lens equation. This matrix may be written as 



60 



(2) 



■K — Jl 
-72 



—72 
K + 71 



where tp(6) is the lensing potential, a two-dimensional projection 
of the gravitational potential along the line of sight (e.g. lSchneiderl 
120061) . These equations define the two components of weak shear 
71 and 72, and the convergence «, which is a measure of the pro- 
jected matter density. 

However, the assumption of gradual variation in lens proper- 
ties across the sky is not always justified, especially in very dense 
regions or those with significant dark matter substructure along the 
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line of sight. In these situations the lens transformation is more ac- 
curately described by 



A6>- ~ AijAOj 



—DijkA6jA9k, 



which is simply the expansion of {T} to second order, where 

d 2 e[ 



D ijk = 



89 ,89k 



= d k Aij = -didjd k ijj. 



(3) 



(4) 



Equation <[3j describes the lensing distortions known as flex- 
ion, which skew galaxy light distributions and lead to weak lensing 
arcs, and which may be described upon an image 1(0) using the 
conservation of surface brightness under lensing 

1(0) = J (s) (6»') =&\e' c + AO'). (5) 

Employing equation ^ in this surface brightness transformation 
takes weak lensing one order closer to the fully generalized non- 
linear treatment. By expan ding the surface brightnes s as a Taylor 
series and substituting {3j, iGoldberg & Baconl d2005l) showed that 
the lensed surface brightness of a galaxy may be approximated as 



7(0) ~ J (s) (6>)+ 



(A - I)ijA0j + -D %jk A9jA6 k 



dj (a) (9),(6) 



an expression which is useful in the construction of weak shear 
and flexion estimators usi ng shapelets dRefregierl2003l ; lBacon et all 
120061: IMassev et alj|2007t) . 

For its clarity and convenience we will often employ the com- 
plex notation introduced in B06. The complex gradient operator on 
the sky plane is defined as: 



8 = di +182- 



(7) 



It is shown in B06 that using this notation the convergence and 
complex shear 7 = 71 + i72 may be written as: 

1. 

* = 2 ( 
1 



7 



-d*dtp, 

-ddip = He 211 *, 



(8) 
(9) 



which neatly encapsulates the spin-2 rotational symmetry proper- 
ties of the shear. Taking a further complex gradient, we may define 
two more complex fields: 



T = -d*ddtp = l^le 1 
G = ^8884- = \Q 



3i0 



(10) 

(11) 



referred to as the first flexion (spin-1) and second flexion (spin- 
3) respectively. Using equation 10 we may write Dij k entirely in 
terms of the components of T and Q, as follows: 



D 



ijX 



D 



ij2 



3Fi + Q x 

T2 + G2 
T2 + G2 

Tx - Gx 



T2 + G2 
Tx~Gx 

Gi 
G> 



(12) 



Tx 
3J" 2 



Equations (|2j, © and j 1 2b allow the formulation of practical es- 
timators for 7, T and G- In the following Section we outline how 
shapelets may be used to make estimators of shear and flexion. 



3 SHAPELET MODELLING 
3.1 Shapelet basis sets 

The underlying c once pt of the shapelet approac h, as introduced by 
iRefregierl l l2003h and lBernstein & Jarvisl d2002l) . is the expression 



of an object's surface brightness as a sum of orthonormal, two- 
dimensional basis functions: 



/w= E E 



. B„ 



(13) 



n\— n<i— 



The choice of basis functions is free in general, but the Cartesian 
shapelet basis set is defined by the basis function 



H ni (9 1 //3)H n2 (9 2 /l3)e-^ 2 ^ 2 
2(™i«2) / g^/ 7rni7l2 



(14) 



where H n ./ X ) is a Hermite polynomial of order m, and the impor- 
tant free quantity (3 is the scale size of the shapelet basis set. We 
refer to the sum of the two parameters nx and n,2 as the order of the 
shapelet basis function, and will generally truncate shapelet models 
to some limiting order n max such that nx + 712 < n max . 

The formalism of polar shapelets, introduced by 
iMassev & Refregierl |2005j), is closely related to that of Cartesian 
shapelets. Instead of the basis set defined by Equations U3\ and 
dl4l l. polar shapelets express the object surface brightness f(6) as 



f(0) = f(e,4>) = J2 E U,™P™Mo,frP), 



(15) 



n — m— — n 



where 9 is the modulus of the complex sky position vector 9± +i#2, 
and 4> = arctan (62 /Ox)- The polar shapelet basis function s, which 
we label P„. m , are defined bv lMassev & Refregierl d2005h as 



/3\™\- 



(n-|m|)/2 \ p2 ' 



7T[(n-|m|)/2]! 

Q2 



-eV 2 / 3 -lm0 



using the follo wing definition of the a ssociated Laguerre polyno- 
mials (see, e.g jArfken & Weberll2005h : 



L q p (x) = 



x~ q e x d p 



P 



dxP 



(17) 



Each separate member of the basis set is uniquely described using 
the two integers n and m, with n > and |m| < n. 

Both the Cartesian and polar sha pelet basis set ha ve relatively 
simple b ehaviour under convolutio n dRefre gier 2003) and decon- 
volution (Refregier & Baco r] |2003h : this makes them particularly 
suited for correcting images for the effects of an instrumental point 
spread function (PSF). We now describe the methodology required 
for this correction. 



3.2 Image deconvolution using shapelets 

Within the shapelet framework, there are two possible methods 
with which to correct galaxy images for the effects of the op- 
tical system. Both approaches begin with the construction of a 
shapelet model of the PSF g(0) at the location of each galaxy: 
this model should be as accurate as possible and will include any 
variation of the PSF across the imag e plane of the instrument (e.g. 
Jarvis & Jair]|2004l: iHoekstral 120041 : iRowel |2010| : iHevmans etaU 
20121 : iKitching et alJl20l3T The model may also need to include 
some treatment of time dependent effects (see, e.g . jHevmans et al.l 



120051 : ISchrabback et alj|2007l : | Rhod es et al 1120071) 

The deconvolution metho d used in this work is that proposed 
by IMassev & Refregierl J2005t) . which is implemented within the 
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shapelet software package made available by these author^ in 
this approach, the deconvolved shapelet coefficients f ni ,n 2 are es- 
timated by "forward" convolving the shapelet basis functions with 
the PSF model in advance, creating a new basis set which we label 



D ni , n9 (0-,l3)=g(0)*B ni , na (0;p), 



(18) 



with an equivalent expression for the case of the polar shapelet ba- 
sis functions P„. m (0; /?). Fitting the data h(9) with this new basis 
set D nit n 2 , one returns a deconvolved shapelet model as follows: 



h{0) = g{0)*f{9) 
= 9{0) 



(19) 



oo oo 



_ri\ =0 n-2=0 



oo oo 



= E E fni,n 2 \g(O)*B ni , n2 (0;P)) 

77,1= 71-2—0 



OO oo 



— E E /"i,«2- D »i,n 2 (^;^), 

ni — T12 — 

As can be seen by comparison with Equation dl3t , the returned 
shapelet coefficients f ni jU2 will reconstruct the deconvolved image 
when they are used with the original basis set B ni ,„ 2 (9; j3). 

There are obvious caveats to this approach, particularly that 
the convolved basis set D ni , n2 (8; j3) will in general no longer 
be orthogonal. However, errors due to this fact are small wher- 
ever the scale size of the galaxy image is larger than that of the 
PSF jMassev & Refregierl 120051) . The alternative shapelet decon- 
volution approach is that descri bed by Refregier & Bacor] J2003h 
and developed in some detail by iMelchior et alj d2009h " involving 
the inversion of a 'PSF matrix' that describes the transformation 
between the shapelet model coefficients of / and h. This PSF ma- 
trix is large and may be sparse, even desp ite efficient truncation 
jRefregieil 12051; iRefregier & BacorJ [20031) . causing the lattei au- 

~~m 



thors and 



to argue agai nst its inversion 



l lMassev & Refregiei 12005 
as a slo w and potentially unstable process. However. IMelchior et al .1 
J2009h find this not to be necessarily the case for sufficiently small 
and lightly-structured PSF models, and so argue in favour of a mod- 
ified inversion scheme for computational efficiency. 

For ACS data we find (see Section l4~3l > that a simple, low n ma x 
model is insufficient to describe the PSF. This makes matrix inver- 
sion schemes problematic, and thus for this .HST-focused analysis 
we will test a shape measurement pipeline based upon the soft- 
ware o f lMassev & Refregiej $200$) . using the forward deconvolu- 
tion method described above. 

The estimation of deconvolved galaxy images using shapelet 
modelling is important at two stages in this Paper. In Section [6] 
shapelet models will provide estimates of the weak lensing sig- 
nal in samples of simulated images. First however, the technique 
is used to estimate the deconvolved shapes of a sample of galaxies 
from which simulated images will be created. Shapelet models are 
convenient for creating such simulations: they are trivially rotated 
and inverted, and as they provide an estimate for the full surface 
brightness profile of a galaxy they allow a full range of distortions 
(including flexion) to be simulated. In the following Section we 
describe the construction of a shapelet galaxy weak lensing simu- 
lation using a set of deep ACS galaxy images. 



http://www.astro.caltech.edu/~rjm/shapelets/ 



Table 1. SEXTRACTOR configuration parameters used to detect galaxy 
objects in both the HUDF (Section [4} and in the simulated ACS images 
(Section[6); these values a re the same as use d fo r the GEMS survey gal axies 
in the two-pass strategy of|.Rix et al] j20()4 and lCaldwell et alj 1200 8|) . 



Configuration Parameter 


Cold Sample 


Hot Sample 


DETECT.THRESH 


2.30 


1.4 


DETECTJMINAREA 


100 


37 


DEBLENDJVlINCONT 


0.065 


0.060 


DEBLENDJMTHRESH 


64 


32 


BACK.SIZE 


214 


214 


BACK_FILTERSIZE 


5 


5 



4 SHAPELET MODELS OF THE HUBBLE ULTRA DEEP 
FIELD 

In order to simulate lensing data, we require a set of real galaxy 
images which are used in turn to create a set of galaxy shapelet 
models which we will refer to as the 'starter set'. The real images 
are obtaine d from the publicly-av ailable Hubble Ultra Deep Field 
(HUDF: see lBeckwith e t al. 2006 for a detailed description), a mul- 
ticolour galaxy survey image composed of a single ACS pointing, 
with a 10 6 s total integration time over four broad spectral bands: 
F435W (-B435), V606W (Veoe), F775W (1775) and F850LP (z 85 o)- 
The highest redshift objects are only visible in the 1775 and zsso 
bands, but the Veoe filt er provides good sens itivity to objects at 
redshifts lower than ~ 4 dBeckwith et al.ll2006l) . As this is also the 
ACS filter used in the GEMS and STAGES lensing studies (du e to 
its rich source number density, see iHevmans et all 20051 [20o3) . we 
choose the Veoe filter image for the construction of the shapelet 
starter set, as representative of typical lensing source galaxies. 

We model the HUDF galaxies using version 2.2 of the 
shapel ets software package, presented by Massev & Refregier! 

and described in practical detail The mod- 

elling of these galaxies, including PSF deconvolution, is a multi- 
stage process: catalogue creation, postage stamp image creation 
and PSF modelling all precede the construction of the deconvolved 
shapelet catalogue for the starter set. These processes are now de- 
scribed in turn. 



4.1 Star and galaxy selection 

The starting point for building accurate models of galaxy images is 
an obj ect catalogue. A primary science goal of the lBeckwith et al] 
( 120061) analysis was a catalogue of full, multicolour photome- 
try in the HUDF, but we require only reliable object detection 
in t he Veoe in order to build galaxy models in the same filter. 
The lMassev & Re fregier ( 2005|) shapelet software package requires 
certain input parameters fo r each object t h at we re not all pro- 
vided in the catalogues of iBeckwith et alj l l2006l) . We chose to 
construct our own catal ogue from just the V6O 6 data, using the 
SEXTRACTOR software teertin & Arnouts|[l9961 : version 2.5.2) in 
single-image mode. Using this software successfully requires a 
number of choices regarding configuration parameters, and we now 
describe the strategy adopted for object detection in the HUDF Veoe 
data. We note that this same strategy will also be used in Section l6~T1 
as part of the lensing measurement pipeline being tested, but for 
simulated imaging data at a much shallower depth when compared 
to the HUDF data used to provide galaxy models. 

We adopt a two-pass SEXTRACTOR debl ending strategy 
when constructing a source catalogue (see, e.g.. iRix et alj|2004t 
Leauthaud et al. I l2007l : ICaldwelT et al.l [2008h . Two catalogues are 
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Figure 1. Size-magnitude diagram for SEEXTRACTOR-selected objects in 
the HUDF Vfj06 science image, showing the stellar locus and 30 selected 
stars (star-shaped points). The dashed lines show the size and magnitude 
limits used to define the galaxy sample chosen for shapelet modelling for 
inclusion in the simulation starter set. 



created, one with a low detection threshold so as to pick out as 
many faint objects as possible, and one using a more conservative 
detection strategy so as to limit the over-deblending of bright ob- 
jects; these catalogues will be referred to as 'hot' and 'cold', respec- 
tively. These are created using exactly the same input parameter 
values as used bv lCaldwell et al.1 d2008l) to create the hot and cold 
samples of objects in the GEMS Veoe science tiles, summarised in 
Table □ 

All cold detections are then combined with non-overlapping 
objects in the hot catalogue. We define a hot object as overlapping 
if its centroid lies within an ellipse of semi-major axis 5.5 x a and 
semi-minor axis 5.5 x b, where a and b are the SEXTRACTOR- 
output semi-major and semi-minor axes, respectively, and this 
larger ellipse is aligned with that defined by SEXTRACTOR. This 
factor of 5.5 was found, by visual inspection of the HUDF and 
segmentation maps output by SEXTRACTOR, to provide a suitable 
compromise between cold object overdeblending and the erroneous 
removal of hot objects. Finally, a mask is applied so as to exclude 
detections from the boundary regions of the CCD image. The com- 
bined hot/cold catalogue then contains a total of 8203 objects, cor- 
responding to ~900 detections arcmin~ 2 . 

To select stars we use the fact that the FLUX_RADIUS param- 
eter found by SEEXTRACTOR is typically constant, irrespective of 
flux, for stellar images. This radial profile of the PSF will typically 
vary little across a given single ACS tile (the same is therefore also 
true of the object Full- Width at Half Maximum: FWHM). This al- 
lows the stars to be easily identified via a straight locus on a size- 
magnitude diagram, such as that showing all the masked HUDF- 
selected objects in Figure [T] The locus chosen in the figure gives 
a total of 30 stellar images from which to build a model of the 
HUDF PSF, avoiding the confused region at greater magnitude and 
saturated images at lower magnitude. This low number of stars is 
to be expected: the HUDF was specifically chosen as a direction 



out of the plane of the Milky Way containing as few stars as possi- 
ble. These 30 stars will be decomposed into shapelet models of the 
HUDF PSF in Section|43l 

To isolate galaxies for shapelet modelling in Section l4l4l those 
to be included in the simulation starter set, we cut for objects with 
20.0 < V 6 06 < 29.0 and 2.2 < FLUX.RADIUS < 150.0 (in pix- 
els). These cuts can also be seen in Figure Q] (barring the large- 
radius cut at FLUX_RADIUS = 150.0). There are then 4128 galaxies 
from the original masked sample that make these cuts, correspond- 
ing to approximately 460 galaxies arcmin~ 2 . It is this extremely 
deep sample that will be used for generating simulated galaxy im- 
ages later in the analysis, although of course many of these ob- 
jects will be lost in noise when simulating shallower data than the 
HUDF. 



4.2 Postage stamp image extraction 

The decomposition of stars and galaxy images into shapelet models 
must be preceded by the creation of 'postage stamp' images of each 
object in the catalogue, cropped around the object in question, and 
masked for neighbours. A postage stamp is also made containing a 
map of the noise and sky background in the same vicinity. 

For each object, the shapelet software draws circular postage 
stamps centred on the SEXTRACTOR-measured centroid. The ra- 
dius of this circular image is the integer number of pixels closest 
to a value rps, defined in terms of the user-specified shapelet in- 
put parameter NFWHM as rps — NFWHM x a + 4, where a is the 
SEXTRACTOR-output semi-major axis of the object. The name of 
the parameter NFWHM appears to be from an earlier incarnation of 
the software in which rps was defined in terms of the FWHM. We 
choose NFWHM = 6 in this analysis, the default value of 5 being 
found to be often insufficiently large to allow the shapelet mod- 
elling of the extended outer profiles of deep, space-based images. 

The default shapelet software flags as a modelling failure any 
objects for which the outer boundary of the shapelet model closely 
approaches the edge of the postage stamp. It was found that while 
the drawing of overly large postage stamps (NFWHM > 7) around 
every object was computationally prohibitive at the shapelet mod- 
elling stage, smaller postage stamps led to an unacceptable number 
of model failures due to the extension of light profiles beyond the 
postage stamp boundaries. An algorithm for iteratively redrawing 
the postage stamp in the event of such model failure provided an 
efficient solution to this problem, and is now part of the shapelet 
software, documented and available for download online at the web 
location given in Section [3] In the iterative prescription used for 
this analysis, model failures due to postage stamp outgrowth are 
resubmitted using a new postage stamp that is increased in size 
by a factor REDRAW .FACTOR = 1.3. This process is repeated 
up to a maximum of MAX_N_REDRAWS = 5 times, after which 
a catastrophic failure is flagged. In tests, these parameter choices 
were found to give a better compromise between modelling suc- 
cess rates and computation time than other values tried: the number 
of floating point operations required to generate images approxi- 
mately varies oc r| s , so REDRAW FACTOR and MAX JM .REDRAWS 
cannot be made arbitrarily large without computational cost. 

After the drawing of a postage stamp around each object of in- 
terest, a fundamental consideration in the shapelet approach is then 
the masking of other, nearby objects: failure to do this well will of- 
ten result in the partial modelling of nearby objects as part of the 
object of interest. Unlike in the KSB approach, image pixels at a 
distance from the object centroid are not explicitly down-weighted 
increasing the importance of careful masking. The shapelet soft- 
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Table 2. Shapelets software input parameters. Upper section: input param- 
eter settings that influence the construction of postage stamps. Lower sec- 
tion: parameters that control the shapelet modelling and deconvolution. No 
default value is indicated for parameters added as a modification to the pub- 
licly available software. 



Shapelets Input Parameter 


Chosen Value 


Default Value 


NFWHM 


6 


5 


REDRAW .FACTOR 


1.3 




MAX_N -REDRAWS 


5 




MASK_NEIGH 


4.0 


2.75 


NEIGHBOUR 


lt.O* 





SKY 


1 





N_MIN (minimum n max ) 





2 


N_MAX (maximum n max ) 


20 


20 


THETA_MIN_GEOM 


0.5 1 ", 1.0* 


0.2 



t values chosen for modelling HUDF galaxies 

* values chosen for modelling the ACS simulated galaxies 

ware therefore constructs over each object an elliptical mask de- 
fined with semi-major axis MASK_NEIGH x a and semi-minor axis 
MASK.NEIGH x b, aligned with the SEXTRACTOR-defined object 
ellipse. The factor MASK.NEIGH = 4 was found to give better re- 
sults than the shapelet default value of 2.75: in the default setting 
there were often portions of the outer galaxy light profile that were 
unmasked and clearly visible. These caused an enhanced rate of 
catastrophic modelling failure for the central galaxy of interest in 
such postage stamps. The larger value MASK.NEIGH = 4 provided 
a significant reduction in such failures without generating unac- 
ceptable numbers of cases where the nearby mask obscures the ob- 
ject of interest in the postage stamp. 

An inverse- variance noise weight map and an estimate of the 
sky background level are then made via analysis of blank sky pix- 
els in the postage stamp, i.e. those unmasked by the central object 
or a neighbour. The sky background can be subtracted from the 
image postage stam p by fitting a choice of surfaces to blank sky 
pixels jBergej|2006h . For the HUDF Veoe image only a very small 
amount of residual sky background variation was found and the re- 
moval of a simple constant sky level from each postage stamp was 
sufficient: this is achieved by setting the shapelet input parameter 
SKY = 1. In order to make the noise map, the root mean square 
(RMS) blank sky pixel value is estimated and then used to create 
a constant, inverse- variance weight postage stamp of the same di- 
mensions as the science image. 

There is a further choice in how noise values are assigned 
for pixels corresponding to masked neighbours, a choice con- 
trolled by the input parameter NEIGHBOUR. For the default value 
NEIGHBOUR = these pixels are assigned zero values in the 
inverse-variance weight map, and are therefore not considered in 
the shapelet modelling. For NEIGHBOUR = 1, the pixels are as- 
signed the same weight as elsewhere in the noise map and set to the 
background level in the science image. Although it is arguably bet- 
ter not to include these masked pixels in any fit for shape inference 
purposes, for the purposes of building a simulation galaxy starter 
set from HUDF images it was found that NEIGHBOUR = some- 
times caused large negative flux patches in shapelet models for un- 
constrained regions of the image. So as to build as physically repre- 
sentative a starter set as possible we therefore set NEIGHBOUR = 1 
for the modelling of HUDF objects, but retain NEIGHBOUR = 
when later testing the shapelet lensing analysis on the derived, nois- 
ier, simulated galaxies. 



UDF average PSF 




-0.2 -0.1 0.0 0.1 0.2 

x [arcsec] 

Figure 2. PSF pattern created from the weighted average of 30 shapelet 
models of selected stars in the HUDF Vgoe science image. The greyscale is 
linear in surface brightness whereas the contours are logarithmic. 

Finally, after these choices for the construction of the postage 
stamps, the trimmed, masked, sky-subtracted science images and 
inverse-variance noise weight maps are then ready to be supplied 
to the shapelet decomposition and modelling routines as described 
in iMassev & Refregien d2005h . We now describe the use of this 
software to model first the stars in the HUDF, and then the de- 
convolved, high-resolution galaxy images that will be used in the 
flexion and shear simulations. 

4.3 Modelling the HUDF PSF 

Modelling the PSF in the HUDF is important. A good PSF model 
allows an approximate deconvolution of the PSF from galaxy im- 
ages, desirable for generating a starter set of models which ac- 
cura tely reflects real galaxy p roperties for space-based data (see, 
e.g.. iMandelbaum et alj|2012l who demonstrate a novel approach 
to simulating convolution-corrected galaxy images). An alternative 
approach is simply to model the ACS YSF -convolved galaxy im- 
ages of the HUDF and use these as the starter set: such an approach 
would be acceptable if the final use of the starter set was in simulat- 
ing observations with far larger PSF sizes, su ch as for ground-base d 
data (this was the approach taken in STEP2: lMassev et alj|2007bh . 
However, when simulating ACS data such an approach would lead 
to an unrealistic size distribution for the small, fainter objects that 
are the most important carriers of weak lensing information. Fainter 
galaxies, being not much larger than the ~ 0.1 arcsec PSF typical 
in ACS images, would be noticeably over-smoothed, too large, and 
with an unrepresentative radial profile. We therefore correct for this 
image blurring as much as possible. 

Taking the 30 stellar objects selected as described in Sec- 
tion |4.1| (see also FigureQJ we first create masked image and noise 
postage stamps as described in Section 14.21 We then construct 
shapelet models of each star, choosing fixed values of /3 = 1.80 
and n max = 20. The fact that these values are fixed, and not al- 
lowed to vary as under the amoeba-driven optimization described 
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by|M assey & Refregierl < t2005h is important: we wish to combine 
shapelet models to create an average PSF for the field, so a fixed 
n m ax and /3 are practical as they allow simple linear coaddition of 
models. The shapelet software also outputs the diagonal entries of 
the covariance matrix for all shapelet coefficients, and so these are 
used to combined the models for each star in an inverse-variance 
weighted average, giving the resulting PSF model seen in Figure[2] 
That this model incorporates no spatial variation across the 
HUDF field-of-view is not a great concern for the purposes of cre- 
ating a realistic starter set for lensing simulations. Although the 
bulk effect of image convolution will be largely corrected for, it 
does mean that some residual ellipticity and flexion will remain in 
the starter set images due to the residual anisotropy variation in 
the HUDF PSF. These faint distortions will therefore be retained 
in the starter galaxy models, but as will be seen in Section[5]these 
starter models are then randomly rotated, inverted and further dis- 
torted before being used in simulations. Within these simulations 
they are then lensed, reconvolved with a new PSF and significantly 
noise-degraded. Distortions due to variation in the original HUDF 
PSF will not have a significant, coherent impact on the final galaxy 
images at the level of measurement possible for simulations of this 
size. 

4.4 Modelling the HUDF galaxies 

In constructing deconvolved (i.e. corrected f or PSF convolution) 
shapel et models of the HUDF galaxies, the iMassev & Refregierl 
d2005l) software takes as its inputs the object catalogue constructed 
as described in Section |4~T1 the accompanying postage stamp im- 
ages and noise maps, and the shapelet model of the PSF. The best- 
fitting shapelet models output by the code will make up the starter 
set that is used to create simulated galaxies. 

However, as was necessary for the modelling of the HUDF 
PSF, choices must be made for input parameter values that gov- 
ern how these best-fitting shapelet models are selected. The most 
important of these are summarised in Table[2] and now described. 

Unlike when modelling the PSF, for constructing galaxy mod- 
els we allow n max and /3 to vary. These quantities are cho- 
sen by an amoeba-drive n optimization process (described by 
IMassev &~R efregier 2 0051 ; we note that other implementations of 
the shapelet m ethod do not necessarily allow n max to vary, e.g., 
lKuiikenll2006h . However, the shapelet software does require the 
choice of a limiting maximum value of n max = N.MAX, and a 
lower starting point N_MIN from which to begin the amoeba search. 

The choice of N.MAX is largely motivated by computing 
resource constraints. In general, a significant fraction of over- 
all processing time is spent modelling a very small subset of 
large/bright galaxies with complex structure. The time taken and 
memo ry required to model a gi ven object increases roughly as 
7in lax dMassev & Refregieill2005l) . Above n max = 20 the calcu- 
lation of pre-multipliers for the shapelet basis functions (see equa- 
tion [76} causes numerical overflow in the unsigned, 64-bit inte- 
gers the shapelets software uses for the calculation and tabulation 
of factorial terms in the numerator and denominator of P n , m . For 
n max > 20 the shapelets software therefore calculates these num- 
bers as double precision floating point numbers, as and when they 
are required. This leads to significant processing overheads, albeit 
ones which might plausibly be averted with some changes to the 
internal processing design of the shapelet software itself. 

However, the largest, brightest and best resolved galaxies are 
also the most likely to be local and only very weakly lensed in 
real data, and this reduces the motivation for dedicating significant 
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Figure 3. Histogram of reduced \ 2 values for the HUDF deconvolved 
galaxy models; the dashed line is the cut introduced to exclude very poorly- 
modelled galaxies. 



amounts of development effort and processing time to producing 
high-order shapelet models of this population. We therefore choose 
N_MAX= 20, which nonetheless allows a high degree of galaxy 
substructural complexity to be realised (e.g. Figure[7}. We choose 
N_MIN = so as to give the shapelet code freedom to model galax- 
ies with a simple circular Gaussian if there is no strong statistical 
support for a more complex model. 

The other input parameter chosen to differ from the default 
choice was THETA_MIN_GEOM = 0.5. This parameter s e ts a lo wer 
limit on the quantity defined by IMassev & Refregierl d2005l) as 
#min = /?/\/iimax + 1, which is the minimum geometric scale on 
which the shapelet model varies. The requirement upon the shapelet 
modelling a moeba of only ac c epting models with 9 m in > 0.5, also 
suggested by iMelchior et alj J2007I) . prevents excessive sub-pixel 
modelling that can lead to unphysical 'ringing' in shapelet models 
of galaxy images. 

Given these chosen input parameters, the postage stamps, and 
the model of the HUDF PSF, the shapelets software was then used 
to find a best-fitting shapelet model of each HUDF galaxy. A total 
of 3867 of the 4128 input HUDF galaxies (94%) were successfully 
modelled. 

On examination of these 6% of catastrophic modelling fail- 
ures, it was found that they were caused either by: 



(i) Galaxy images lying too close to the edge of the HUDF or 
other objects in the field (being thus automatically rejected and 
flagged in the output shapelet catalogue); or 

(ii) A population at faint magnitudes which appear to suffer 
from either star-galaxy confusion or general modelling degeneracy. 

The great majority of this second type of catastrophic failures were 
for objects with Veoe > 27, which is currently beyond the realm 
of plausible lensing measurement with even deep imaging surveys 
such as GEMS or COSMOS. Without a clear means of further re- 
ducing this failure rate, it is tolerated given the fact that it predomi- 
nantly effects objects beyond the detection limit of the simulations 
(See Section [5~6l >. 

An additional cut was imposed on the sample, based upon the 
value of the reduced \ 2 output for each galaxy model. A histogram 
of this output statistic is shown in Figure[3]for the HUDF galaxies. 
A cut of reduced % 2 < 5, beyond which the chances of a good 
fit are vanishingly small for the high degree-of-freedom shapelet 
models used, removed a further 153 very poorly-modelled galaxies 
not identified by shapelet flags or other indicators of total mod- 
elling failure. The remaining 3714 PSF-corrected galaxy models 
form the galaxy starter set that are now used to construct sim- 
ulated lensing survey images with realistic galaxy morphologies. 
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This corresponds to a high density of approximately 410 galax- 
ies arcmin -2 , giving considerable freedom to realistically simulate 
shallower HST data than the HUDF (an important caveat is the lim- 
ited size of the sample itself, which we discuss in the following 
Section). 



5 SIMULATIONS OF ACS DATA 

5.1 Source galaxy images 

The galaxy models in the starter set described in the previous Sec- 
tion are the fundamental data used to generate our simulations. 
However, there is additional processing required before models 
may be used to accurately simulate flexion and shear measurement. 

A first problem is that the starter set itself contains the sig- 
nature of gravitational shear and flexion due to the matter structure 
along the line of sight of the HUDF, along with any residual, uncor- 
rected HUDF PSF anisotropy. These signals are largely eradicated 
by the random rotation and inversion of galaxy models. These two 
transformations can be performed by simple analytic manipulations 
of the galaxy shapelet mod els in the polar shapelets formalism of 
IMassev & RefregiejJiool) . and so this process is fast and accurate. 

A second problem is that the starter galaxy set represents a 
limited sample of galaxy morphologies. This may be alleviated 
to an extent by introducing small random perturbations to the 
shapelet models of the galaxy ima ges (e.g., |M assey et al. | |2004 
iMeneghettietal .1120081) . We follow Meneghet ti et al .1 d2008t) and 
add a Gaussian random variable N(0, a ni ,n 2 ) to the real and 
imaginary part independently each shapelet coefficient fm,n 2 - We 
choose <J ni ,n 2 — 0.1 x |/ni,n 2 |> which we found ensured that 
galaxy shapes are randomized in a way that did not introduce a 
high inci dence of noticeably unphysical features such as negative 
flux (see IMassev et al"]|2004l for a discussion of over-randomized 
galaxy models). 

Randomly rotated, inverted and perturbed galaxy models are 
then a suitable population of source galaxies, and can be assigned 
to locations in the final output images. U n like the GREAT08 o r 
GREAT 10 simulations dBridle et al.ll2010l : iKitching et alj|2012h . 
we choose to assign galaxy images to random locations upon the 
output science tile, rather than on a regular grid. Although this 
complicates the analysis it has the benefit that the impact of object 
detection and deblending is included in our pipeline tests. These 
effects are important as they likely contribute to noise on shear and 
flexion estimates, one of the primary topics of investigation in this 
work, enhancing the value of any noise forecasts based on the re- 
sults. 

5.2 Applying shear and flexion: shapelet transformations 
versus raytracing 

In the previous Section we arrived at a set of shapelet models suit- 
able for use as source galaxies. The desired weak lensing distor- 
tion, a coordinate transformation governed by the lens equation, 
must now be applied to these galaxies. There are two approximate 
means by which this coordinate transformation may be applied in 
practice: 

(i) In 'shapelet-space': shear and flexion distortions may 
be represented by an infinite sum of transformations upon the 
shape l et coefficient values in the galaxy model itself dRefregieil 
120031 : IMassev & Refre gier 2005|; IMassev et all |2007|) . Applying 
successive shapelet transformations therefore allows the lensing 



distortion to be represented with arbitrary precision. The lensed 
shapelet model may t hen be analytically integrat ed across square 
pixels as described bv lMassev & Refregien J20051) . 

(ii) In 'real-space' : the shear and flexion distortions may be rep- 
resented directly using the lens equation 1(0) = /' s) (0') to ray- 
trace from the regular grid of image pixels back to an irregular sam- 
pling of points on the source plane. This is possible as each shapelet 
model encodes the surface brightness at all points in the source 
plane. 

Previ ous shapelet-based simulation methods have p referred the for- 
mer dMassev et al]|2007bl : iMeneghetti et al.ll2008l) for weak lens- 
ing transformations. For our simulations, designed to mimic space- 
based flexion and shear observations as closely as possible, we 
choose the latter real-space method for reasons that we now de- 
scribe. 

In the shapelet-space approach it is necessary to apply multi- 
ple shapelet transformations to each galaxy model, each operation 
capturing successive terms in a series expansion representation of 
the true distortion, to ensure that the final distortion approximates 
a real shear or flexion accurately. Doing this is important if subse- 
quently we wish to test sh apelet measurement metho ds fairly. As 
shown by equation (41) of dMassev & Refregiej2005t) . performing 
a shapelet shear transformation accurate to first order in 7 increases 
the order n max of the model from n max to n max + 2, resulting in 
2n max + 5 extra coefficients overall. At third order in 7, the min- 
imum that should be considered for precision work where the ap- 
plied shear 7 may approach 0. 1, the order of the model is increased 
to n ma x + 6, an extra 6n max + 27 coefficients in total. 

The situation is worse for flexion. iMassev et al. (2007) 
showed that a first-order approximation to gravitational flexion in- 
creases the order of a shapelet model to ?i max + 3. Performing a 
shapelet flexion transformation that is accurate to third order there- 
fore increases n max by 9, meaning that a modest n max = 12 
shapelet galaxy requires a far more complex n max = 21 model 
once flexed. The shapelet software suffers extreme memory de- 
mands and reductions in speed beyond n max = 20, and so this 
method of introducing accurate distortions becomes prohibitively 
slow. Unfortunately, re-truncating the shapelet model back to a 
more manageable n max is not a simple solution to these problems. 
This degrades the congruence between the exact and shapelet- 
approximated lensing transformations in a way that is difficult to 
describe, as it varies with both the distortion applied and the sur- 
face brightness distribution of each individual object. 

Also of concern is the action of the lensing transforma- 
tions upon the j3 scale parameter. Lensing distortions magnify im- 
ages, which for the linear-order shear and convergence transfor- 
mations can be simply represented by a small change in /?. With 
flexion a fundamental complication with ft arises: as shown by 
Schneider & 13 d2008l) . the determinant of the lacobian matrix 
varies as a function of position for flexion-order lensing. This vary- 
ing transformation in the area element cannot be reproduced via a 
single rescaling of the /? scale parameter. 

This potentially increases the difficulty of exactly reproduc- 
ing non-linear flexion transformations using the shapelet formal- 
ism. When constructing simulations for flexion measurement we 
must be extremely careful that we are accurately describing the 
distortion if we wish to construct a fair test of current and future 
methods. The real-space, raytracing option listed above is there- 
fore adopted for applying both flexion and shear in the simulations. 
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We now test the level of accuracy that may be achieved with this 
method. 



5.3 Testing shear and flexion raytracing 

When performing real-space raytracing as described in Section [5~2l 
we use the lens equation to map from the image pixel positions 
back to a non-regular grid of positions in the source plane. At these 
points we sample the surface brightness from the shapelet model 
of the source galaxy I (s) . I gnoring the overall deflection of galaxy 
images (equivalent to setting C — 0' c — 0), and using equation 
l[3}, we may write 



I(0) = I {s >(0') = I (s) (A 



'3 + Di ]k 9j9 k /2). 



(20) 



In practical terms the raytracing scheme is simple: we assign 1(0) 
for each desired image pixel position by taking the shapelet 
model value of / (s) at the position 0' = 0'(O, 7, T, Q). 

Square pixels in the image plane do not map to square pix- 
els in the source plane, and when flexion is included with shear 
the mapped pixel boundaries become curved in the source plane. 
This means that it i s no longer possible to use the results of 
iMassev & Refregiej J2005h to perform the exact flux integral of 
1(0) across each pixel. In the simplest approximation one can 
adopt the pixel-centre surface brightness to estimate flux but more 
accurate results can be achieved by upsampling, i.e., by creating 
a higher-resolution image and then summing values at the high- 
resolution subpixel locations to estimate the true integral. This al- 
lows approximation of the exact shear and flexion transformation 
to an accuracy that depends on the degree of upsampling adopted. 

We can only upsample by a finite, preferably integral, factor. 
In order to understand what factor is necessary in our simulations 
we now test its impact upon fundamental lensing measures (e.g. 
image moments) for noise-free ACS galaxy models. 

We first extracted a random sample of 1000 galaxies from the 
HUDF starter set. For each galaxy, a control postage stamp image 
was constructed by integrating across ACS imaging survey-size, 
post-drizzling pixels (0.03 arcsec, as used in the H UDF, COS MOS, 



GEM S and STAGES final s ci ence images; see Beckwi th et al 



2004 iLeauthaud et al.l 120071: ICaldwell et all 120081 : iGrav et al 



200S ) using the exact Massev & RefregierT i 20051) analytic results 



We refer to this 'true' pixelized control image as I^K Simple, un- 
weighted moments for such images may be defined as follows 



S = 



Hi 

(J i j k 



d 2 ei(o), 

d 2 eAe i A0 j I(0), 
d 2 0A0 i A0 j A0 k I(0), 



q ljk i = g I d 2 6A8 i A9 3 Ae k A9 l I(0), 



(21) 
(22) 
(23) 
(24) 



where AO — — C as before. In the noise-free case the un- 
weighted moments can be used to construct complex polariza- 
tion measur e s that provide estimators of shear and flexion (e.g. 
Kaiser et all 1 19951 ; IBartelmann & Schneider! 1200 ll : lOkura et all 
20071) . 



We construct the unweighted polarization 

gn — 922 + 2igi2 




1 2 3 4 5 6 7 

Upsampling ratio r 

Figure 4. Testing the impact of upsampling ratio for generating images 
of lensed (by shear and flexion) galaxies using multiple 'ray-traced' im- 
age samples. Plot of median fractional error in simple unweighted e (top 
panel), / (centre panel) and g (bottom panel) estimators, for a sample of 
1000 galaxies randomly selected from the starter set, with increasing up- 
sampling ratios r as described in Section |531 Each Sfi etc. is cal- 
culated by comparing the upsampling-derived value relative to the exact 
(analytically integrated) value. Th e wide solid er ror bars on each point give 
the standard error on the median iLupton|[T993h . whereas the dashed error 
bars illustrate the typical range of the effect as described by the normalized 
median absolute deviation (NMAD) of the fractional error. 



(see, e.g.. iKaiser et alj[l995 : IBartelmann & Schneidej[200ll) and, 
following lOkura et alj ( 120070 ■ we define 



/ 



9 



h + i/a 



91 + 132 



<7m + <7i22 + i(gri2 + (7222) 

gnu + 1q\\%i + g2222 
gin - 3gi22 + i(3gn2 — 9222) 
gnu + 2gn22 + g2222 



(26) 
(27) 



ei + ie2 



q\\ + g22 



(25) 



as equivalent measures for T and Q respectively. We require that 
any finite degree of upsampling must cause fractional errors in e, 
/ and g that are significantly smaller than those we expect due to 
noise in the final simulation results. 

The first step in the test is to calculate control values e'*', 
and using the moments as defined in equations d2 lb-(|27|> for 
each control image 1^ from the sample of HUDF starter galaxies. 
We then make a series of raytraced estimates I r of by perform- 
ing a numerical integral of the flux across the upsampled pixels. 
We define the upsampling ratio r as the ratio between the linear 
scale of output pixels and that of the sub-pixels in the upsampled 
image: r — 1 is equivalent to simply taking the central pixel value; 
r = 2 is equivalent to four equally spaced sub-pixels; etc. For each 
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galaxy I r we measure a corresponding e, /, and g, and calculate 
the fractional deviation in each component: 

Se e--e (t) 

^1 = e * e i t (28) 

with equivalent expressions for / and g, for r — 1, . . . , 7. In Fig- 
ure 2] we plot the median of these fractional deviations across the 
sample of galaxies, as a function of r. 

As a measure of the typical range for the error due to finite 
upsampling, in Figure [4] we also plot as dashed-line error bars the 
normalized median absolute deviation (NMAD) of the fractional 
error. The NMAD is a robust measure of the width of a distribution 
of some variable x, defined as 

NMAD (a;) ~ 1.4826 x MAD(i) (29) 
= 1.4826 x Median (\x - Median (x)\) , (30) 

where the MAD is simply the median absolute deviation, and the 
constant scale factor 1.4826 ensures that the NMAD of a Gaus- 
sian distribution is approximately equal to its standard deviation, 
NMAD ~ a. Figure |4] shows the NMAD for Sei/et, Sfi/fu and 
5gi/gi at each value of r, to give an idea of the typical range in 
fractional errors for these quantities due to the use of finite element 
upsampling techniques. It should be noted that median statistical 
measures were used because the property of interest, which in this 
case is specifically the fractional uncertainty on an ellipticity or 
flexion measurement, leads to large outlier values when there are 
small values of et, fi or g, in the denominator. Adopting median 
statistics downweights these noisy outliers to reduce overall statis- 
tical uncertainty on the central tendency. As a check on the results, 
simple means were also used and resulted in a consistent (but rather 
noisier) trend. 

The results of Figure|4]are encouraging and as expected: there 
is clear convergence towards smaller and smaller fractional errors 
with increasing r. Practically, the results imply that for the HUDF 
galaxies and an output pixel scale of 0.03 arcsec we need only over- 
sample galaxies by a factor of two in order to give output images 
that are typically accurate in e, / and g to within a factor ~ 10~' ! 
of the correct value (a targ et value motivated by the results of, e.g., 
lAmara & Refrl gier 2008). Interestingly, Figure[4]shows that good 
accuracy may be achieved for e, / and g, even when r = 1, for 
which the error in the idealized measures is often less than 1% of 
the underlying value. 

To conclude, we have ensured that shear and flexion can be 
added to our simulation galaxies with sufficient accuracy for cur- 
rent purposes. We turn now to the problem of convolution with an 
instrumental PSF 



5.4 Convolution 

Once flexion and shear has been applied to simulation galaxies, the 
next stage is to convolve these correctly-lensed objects with a simu- 
lation kernel that mimics a realistic target PSF. The PSF we choose 
for the simulations in this study is based on ACS observations of 
stars from V606 obse rvations in the GEMS survey dRix et al.ll2004l ; 
ICaldwell et aT]|2008h, specifically the 909 stellar objects selected 
as described bv lHevmans et al.U20051) . Shapelet decompositions of 
these objects are created using n max = 20, fi = 1.80 pixels, and 
stacked to make an inverse-variance weighted average GEMS PSF 
model (the same procedure as described in Section |4~3l for mod- 
elling the HUDF PSF). This model can be seen as the leftmost panel 
in Figure [5] 



As described by|M assev & RefregieJ l l2005t) . any image de- 
scribed using polar shapelets may be easily 'circularized' (i.e. cir- 
cularly averaged) by setting the model coefficients f n ,m = for 
all m 7^ 0. The circularized GEMS PSF generated in this way can 
be seen in the central panel of Figure [5] (imaged at a resolution of 
0.015 arcsec pixel -1 ), and in the rightmost panel (at the final out- 
put resolution of the simulations, 0.03 arcsec pixeP 1 ). We choose 
to use the circularized GEMS PSF for the simulations in this Paper 
as its symmetry simplifies the interpretation of lensing measure- 
ment results, while still incorporating a radial profile characteristic 
of space-based data such as that from ACS. Diffraction spikes, such 
as caused by support struts for the secondary mirror, will be lack- 
ing from this adopted PSF model, but the extent to which these are 
successfully characterized by shapelets PSF models is unclear even 
when not circularizing, as in Figure [2] 

Since shapelet models are no longer being used to describe 
galaxies after shear or flexion is applied, the convolution must be 
performed numerically using a pixelized image of this PSF. As a 
shapelet model P SF such as th at in Figure [5] is not formally band- 
limited (see, e.g.. lMarksll2009]) , this therefore requires another em- 
pirical investigation into the effects of finite sampling. 

That this investigation is numerically feasible also illustrates a 
further advantage of not perfor ming the calculat i on usi ng shapelet 
transformations. As shown by iMelchior et al .1 J2009h . an exact 
shapelet treatment of convolution results in a convolved image of 
order n miix given by 

^■max, convolved — ^max, galaxy ~\~ ^max, PSF- (31) 

The PSF in ACS images requires n max , psf > 20, which makes 
describing perfect convolution on even modestly-sized galaxy 
shapelet models extremely expensive. Re-truncating the convolved 
model to the original n max spoils the exactness of the shapelet con- 
volution treatment in a way that varies depending on the shape of 
each galaxy. In contrast convolutions on images, performed using 
the Fast Fourier Transform (FFT, a fast algorithm fo r performing 
discrete Fourier transforms: see, e.g., |Pressetal.ll992h . can be per- 
formed to great accuracy in a fraction of the time. 

To test the accuracy of numerical convolution we use the 
methods developed in Section |531 We again select a random sam- 
ple of 1000 HUDF galaxies from the starter set. For each galaxy 
we construct a set of 7 images, upsampled by a linear factor 
r = 1, . . . , 7, and perform a convolution via FFT using an im- 
age of the circularized GEMS PSF (Figure [5} upsampled by the 
same factor. The convolved image is then summed back to the final 
resolution of 0.03 arcsec pixel , and then each of e, / and g are 
measured. 

Because there is no practical way of generating the 'true' con- 
volved image in this manner as a control (as discussed before the 
exact shapelet treatment is computationally unfeasible), we instead 
take the high upsampling case of r — 11 as the reference point. 
The fractional deviation from this r = 11 is then calculated for 
each galaxy and each r, and in Figure|6]we plot the resulting me- 
dian values and range as in Section l5T3l 

As was found for the shear and flexion upsampling tests, the 
results of Figure [6] are encouraging: an upsampling ratio of only 
r = 2 allows a systematic fractional error in the estimates of e, / 
and g that is typically 0.1% or less. Interestingly, the effect is in 
the opposite direction to that induced when performing raytracing 
lensing transformations. The effect of a finite sampling approxima- 
tion to convolution is to cause a net reduction in these moments of a 
galaxy image, whereas a finite sampling approximation to the lens 
mapping artificially exaggerates these moments. The overall net ef- 
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Figure 5. Left panel: PSF pattern created from the weighted average of 909 shapelet models of selected stars in the GEMS V606 science tiles. Central panel: 
circularized version of this same GEMS PSF made by setting to zero all the polar shapelet coefficients f n ,m with m 0, and shown at the upsampled 
0.015 arcsec pixel -1 resolution used for performing the real-space convolution. Right panel: the circularized GEMS PSF, shown at the final ACS 0.03 arcsec 
pixel -1 resolution for reference. In all panels the greyscale is linear in surface brightness whereas the contours are logarithmic. 




feet of the two approximations will thus, on average, be typically 
smaller than that described by either of Figures[4]or[6]in isolation. 
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Figure 6. Testing the impact of upsampling ratio in numerical convolution 
via discrete Fourier transforms, using FFT. Plot of median fractional error 
in simple unweighted e (top panel), / (centre panel) and g (bottom panel) 
estimators, for a sample of 1000 galaxies randomly selected from the starter 
set, after performing a convolution using discrete Fourier transforms at an 
upsampling ratio r as described in Section l5~4l Each Sej, Sfi etc. is cal- 
culated by comparing the ovesampling-derived value relative to the value 
found for r = 11. Th e wide solid er ror bars on each point give the standard 
error on the median <Lupton|[T993l) . whereas the dashed error bars illus- 
trate the typical range of the effect as described by the normalized median 
absolute deviation (NMAD) of the fractional error. 



Design strategy for distorted, convolved images for the 
ACS simulations 



5.5 



The results of Sections [5,3l & PJ4l now allow informed decisions to 
be made regarding the construction of the ACS simulations to test 
the measurement of flexion and shear with shapelets. Current shear 
measur ement methods are able to measure shear at percent-level ac- 
curacy JBridle et al.ll201u iKitching et al.ll2012h , and it is unlikely 
that flexion measurement will approach this capability within the 
near term. A conservative requirement is therefore that the treat- 
ment of distortion and convolution in the simulations should be 
accurate at the 1CP 3 level in terms of the impact on 7, T and Q 
in simulated galaxies. This also matches stated requirements on 
the estimation of shear fo r an all-sky, space-based survey mission 
dAmara & Refregiej|2008t) . 

We therefore adopt an upsampling ratio of r = 2 when ap- 
plying both the lensing image distortions (shear & flexion) and 
for subsequent convolution with the circularized PSF of Figure [5] 
corresponding to an absolute resolution of 0.015 arcsec pixel -1 . 
It is computationally convenient that both distortion and convolu- 
tion occur at this same resolution. The results of Sections [5.3| & |5.4l 
suggest that errors in the pertinent moments of simulated galaxy 
images will be thus controlled to better than 0.1% on average. 

A remaining issue is the level of input shear and flexion distor- 
tions which should be applied to the simulation images. For shear 
we require successful recovery of the signal due to correlated large 
scale structure, but also for galaxy-galaxy lensing and the stronger 
shears expected around cluster lenses. For flexion, it is unlikely 
that the cosmological signal is measurable in the near future, but 
the galaxy-galaxy signal may be of interest (e.g. B06) and there 
are certainly applications i n the field of cluster reconstruction (e.g., 
iLeonard et al.ll2007l . l201 ll) . 

For the suite of simulations we choose to split our input sig- 
nals into three broad groups exploring a range of distortion signal 
strengths, labelled as 'high', 'mid' and 'low'. These designations 
refer to the magnitude of the input gravitational signal applied, 
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Table 3. Flexion and shear input signal values for the 270 simulated ACS 
tiles, where the first flexion input is given by Fmvut = |-7"inputl e > the shear 
by 7in P ut = |7input|e 2i * and the third flexion input by g input = |Gin P ut|e 3 "* 
(see Section l5~5V 





Set 


|*^""inputl» 1 Ginput 1 


1 Tinput I 


<t> 


Tiles 






arcsec" 1 










1 


0.05 


0.1 


0° 


1-10 


'high' 


2 


0.05 


0.1 


30° 


11-20 




3 


0.05 


0.1 


45° 


21-30 




4 


0.01 


0.05 


0° 


31-40 


'mid' 


5 


0.01 


0.05 


30° 


41-50 




6 


0.01 


0.05 


45° 


51-60 




7 


0.005 


0.01 


0° 


61-70 


'low' 


8 


0.005 


0.01 


30° 


71-80 




9 


0.005 


0.01 


45° 


81-90 



In order to take advantage of the reduction in shape noise that 
can be achieved by combining lensing measurement s from appro- 
priately rotated galaxy images jMassev et ai]|2007bh . we generate 
a further set of 'rotated' simulation images. These are identical to 
those described above except for an additional rotation of 180°, 
90° and 60° (for T, 7 and Q respectively), given to the starter set 
galaxies immediately after step (ii) above. 

This allows the distortion measurements for matched pairs of 
rotated images to be averaged, cancelling the leading order impact 
of the intrinsic galaxy shape upon shear and flexion measurement. 
This also allows the relative impact of galaxy shape noise and im- 
age pixel noise to be compared, which is of particular interest for 
flexion. The total number of simulated ACS tiles in each suite is 
therefore 2 x 90, leading to a full suite of 6 x 90 = 540 simulated 
ACS pointings (~ 100 Gb of data in total). 



chosen to be |7i np ut| = 0.1, 0.05, 0.01 respectively for the shear 
simulations. 

For the flexion simulations it was decided to bring the 'mid' 
and 'lo' samples closer together (to concentrate on exploring sensi- 
tivity to galaxy-galaxy flexion) and extend the high signal some- 
what further to investigate mea surement of values which may 
be found in cluster studies (B06; iLeonard et alj[2007l 1201 ll) . We 
choose I ^"mput I, |<5input| = 0.05, 0.01, 0.005 arcsec -1 for 'high', 
'mid' and 'low', respectively. 

In order to explore any anisotropy in signal recovery (due, 
e.g., to alignment with pixel axes) we split each of the three sets 
into a further three subsets by the angle of orientation of the input 
signal with respect to the image a>axis. Orientations of 4> = 0°, 



30°, 45° were chosen for the input signals T\ 



input 



7m P ut = |7input|e"'' and t/ taput = |Smp U t|e , giving a total of 
3x3 = 9 subsets overall for each of the three lensing distortions. 
These values and choices are summarised in Table [3] We note that 
these values span a predominantly positive range of values in the 
components of 7, T and Q : this asymmetry in the applied signal is 
benign due to the circular symmetry of the adopted simulation PSF 
(Figure[5}. 

Ten survey tiles were then simulated for each image set de- 
scribed above, each with a 30 x 30 arcmin 2 sky coverage area 
(i.e. an ACS pointing) and at an output resolution of 0.03 arcsec 
pixel . These were created by: 

(i) Randomly selecting with replacement from the HUDF starter 
set described in Section|4] 

(ii) Randomly perturbing, rotating and inverting each starter set 
galaxy model as described in Section [5~TI 

(iii) Applying a lensing distortion as prescribed by Table[3]using 
the raytracing method presented in Section [5~2l with an upsampling 
ratio of r = 2. 

(iv) Convolving each distorted image with the circularized 
GEMS PSF of Figure [5] This convolution was performed at the 
image level using FFTs as described in Section \5A\ again using 
r = 2 for the upsampling ratio. 

(v) Placing each simulated galaxy image at a random position in 
the tile. 

In this way 3 x 90 simulated ACS tiles were created for each of 7, 
T and Q. We note that it was decided not to include simulation tiles 
containing both shear and flexion input signals simultaneously. The 
possibility of cross-contamination between the se signals is inter- 
esting (see lViola. Melchior & Bartelma nn 2012), but in this initial 
study we concentrate on examining each signal individually. 



5.6 Correlated noise 

The final ingredient in the creation of simulated images is the 
addition of realistic measurement noise, due to diffuse back- 
ground light and finite photon number counts. An important con- 
sideration for shape measurement is the impact of spatially cor- 
related noise. This is present due to the standard practice of 
combining multiple, dithered ACS exposures to generate sin- 
gle 'science' images using software such as MULTIDRIZZLE 
i lFruchter & Hookl 120021 : iKoekemoer et al.ll2002h . and is generic 
even for more carefully optimized linear combination schemes 
(e.g.. iRowe. Hirata & Rhodesll201 lh . Such dithered images were 
used as the basis for weak lensing measurement in each of GEMS , 
COSMOS and STAGES llRix et :alj [20041; IHevmans et al.l |2005| ; 
Leauthaud et al.ll2007l : IHevmans et alJl2008| - |Caldwell et alj|2008l : 



Gray et al .11 2009) 



We add realistic correlated noise to the simulation tiles de- 
scribed in Section |531 in a novel manner, using a 'noise mosaic' 
image constructed as p art of the imaging reduction of the GEMS 
survey dRix et al.|[2004h . This composite image, which is the size 
of a single ACS pointing, is a mosaic of multiple bl ank sky regions 
in the GEMS Vsoe science images described by ICaldwell et al .1 
J2008I) . Each 0.03 arcsec pixel 1 GEMS science image was gen- 
erated from a dithered combination of three 0.05 arcsec pixel -1 
exposures, each of total duration 2160s. The composite noise im- 
age therefore reflects precisely the noise-correlating effects of the 
GEMS dither and drizzle strategy, which is typical of the reduc- 
tion strategies employed for high-resolution, space-based imaging 
data. The galaxy image magnitude zero points were set to match the 
GEMS data, resulting therefore in simulated survey images with a 
detection threshold around Vaoe ~ 27, matching the GEMS survey 
characteristics quite closely dCaldwell et al.l2008T) . The fact that the 
images include many more simulated objects, buried in the noise 
down to the HUDF starter set cutoff magnitude of Veoe < 29, is 
another realistic feature of the simulations. 

As galaxies are placed at random locations in each successive 
simulation tile, it is possible to re-use this noise image each time. 
Only a small fraction of galaxy images within tiles will consist of 
galaxy models that, by chance, have been repeatedly placed into the 
same location in the ACS pointing-sized noise mosaic (~ 3.53 x 
3.63 arcmin in total area), so the level of unwanted pixel noise 
repetition in the simulation images will be low. In these simulations 
we do not add an additional Poisson noise term to image pixels to 
account for variation in variance due to varying flux. However, as 
the overwhelming majority of simulated galaxies in the catalogues 
are faint relative to the background this can be shown to represent 
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Figure 7. Illustrative 30 arcsec X 22.5 arcsec section of the 86th tile (chosen 
at random) in the simulated, convolved, J-"-lensed images. The greyscale is 
linear in flux. 



is a Moffat profile motivated by ground-based PSFs (and made el- 
liptical for some sets of simulated images), we instead employ a 
circular PSF with a radial profile taken from shapelet fits to GEMS 
stellar images. Deble nding is an issue in our simulations, whereas 
IVelander et alj fcoill) placed galaxies on a regular grid. Finally, our 
simulations employ correlated noise take n from observations in the 
GEMS survey data dCaldwell et al]2008h . 

In Figure [7] we show a small section from one of these tiles, 
illustrating some of the realistic diversity of galaxy morphology de- 
picted in the simulation images. All of these simulation images are 
available by request from the authors. We now turn to a descrip- 
tion of the measurement of shear and flexion in these images using 
shapelets, allowing a calibration of the shear and flexion signal re- 
covery using this method. 



6 TESTING A SHAPELET FLEXION PIPELINE 

We now describe measurements of flexion and shear made from 
the simulations described in Sections [4] & [5] and investigate the 
recovery of flexion and shear as a function of image properties such 
as SNR and galaxy size. 



a very small correction to the images, arguing that the effect can 
be safely neglected in light of other uncertainties in the modelling 
overall. 

One final concern in the use of the GEMS noise mosaic might 
be if there were to be found to be some preferred direction in the 
noise pattern. Visual inspection of the noise map did not suggest 
any such artifact, and it will be seen in Section[6]that there is no ev- 
idence for an additive bias in shear or flexion result s, or a variation 
in multiplicative 'm' biases dHevmans et~ai]|2006t) with polar an- 
gle with respect to the pixel grid. Global anisotropics or preferred 
directions in the noise map might be expected to cause such effects 
in the presence of a circularly-symmetric PSF such as that adopted 
(Figure O, and the fact that we are unable to detect them within 
simulation uncertainties suggests that they exist at a sufficiently low 
level as to not affect the conclusions of this study. 

It should be noted that by reusing the GEMS noise mosaic for 
each tile in this way we are applying the same noise field to both 
galaxy images in the rotated and un-rotated galaxy pairs, as in these 
pairs the simulated galaxies share the same location. However, this 
offers an opportunity to separate the effects of pixel noise in flexion 
measurements from those of shape noise without diluting the for- 
mer. The price paid for this opportunity is a reduction in the overall 
statistical power, by a factor ~ \[2 in the simplest estimate, with 
which measurement bias parameters may be constrained. 

5.7 Summary and comparison to previous flexion simulations 

After the addition of the noise image, the 540 ACS simulation tiles 
are complete. We now briefly summarize the differe nces between 
these simulations and those of IVelander et lb, in case such 

a comparison is of utility to the interested reader. A primary dif- 
ference is that the simulations in this study employ complex mor- 
phologies in the shapelet galaxy models, constructed from HUDF 
imaging data, compar ed to the simpler parametric forms as used in 
JVelanderetal.ll201lh . Our simulations also employ a continuous 
distribution of object sizes and signal-to-noise, taken directly from 
the UDF sample, rather than fixing gala xy sizes and signal-to- noise 
at fixed values of interest. Whereas the lVelander et al j J201 lh PSF 



6.1 Object detection and shapelet decomposition 

To test and calibrate shapelet measurements of shear and flexion, 
we treat the 540 ACS simulation tiles as if they were new tele- 
scope data (with properties such as galaxy positions and shapes 
unknown). The first step is therefore to detect objects in the images 
from peaks in the surface brightness distribution, and we employ 
the same techniques as described in Section |4~T1 to create a cata- 
logue of galaxy objects in the HUDF. We detect galaxy objects in 
the simulations using the SEXTRACTOR software with the parame- 
ter choices given in TableQ] An initial cut of Veo6 < 27 is then ap- 
plied to the catalogues, along with requiring FL UX_RADIUS > 2.4 
and F LUX_AUTO/FLUXERR.AUTO > 10 (see iBertin & Arnoutsl 
1 19961 for descriptions of these SEXTRACTOR outp ut parameters). 
These cuts are motivated by the choices made in iHevmans et al.1 
( l2005h . and result in catalogues containing approximately 64 galax- 
ies arcmin -2 . This figure agrees well with galaxy densities found 
in surveys at similar depth t o that simulated here , such as G EMS, 
COSM OS, & STAGES (e.g. jHevmans et al.l2005l : lLeauthaud etall 
|2007|) . 

Postage stamp images of each detected galaxy object are then 
created as described in Section 14.21 and shapelet models of the 
galaxies are created as described in Section 14.41 except for two 
important differences that we now describe. The first and most 
important difference is in the PSF: for this we use the shapelet 
model of the circularized GEMS PSF (Figure[5j described in Sec- 
tion [5T4] Therefore our calibration of flexion and shear measure- 
ment does not include the potential effects of a poorly-modelled 
PSF, such as might be present working with real data. While 
unrealistic, this simplification will allows any measurement bi- 
ases to be interpreted cleanly rather than being subject to exter- 
nal factors such as poor PSF modelling. The problem of build- 
ing accurate PSF models is important enough to be addressed 
in its ow n right, and this is increasingly reflected in the litera- 
ture (e.g. iHoekstral |2004 Ijarvis & Jainl [20041; Haryis et alj|2008l: 



Paulin-Henriksson et al.l l2008l 120091 : iRowd |20ld : fetching et atl 
201 it Etching et al.. in prep.). 



The second important difference is in the choice of the 
NEIGHBOUR input parameter to the shapelet software, setting 
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NEIGHBOUR = 1. As described in Section ET2l this causes pixels in 
the masked areas of each postage stamp (i.e. those associated with 
nearby objects, or cosmic rays, bad pixels etc. in real data) to be 
given zero weighting at the shapelet modelling stage. Having made 
these changes to the input settings, the modelling provides a cata- 
logue of shapelet coefficient values for each galaxy in the rotated 
and unrotated simulations. 

The galaxies in the rotated and unrotated simulations are then 
matched, treating as pairs all galaxies with centroids separated by 
less than 0.15 arcsec in centroid (five 0.03 arcsec pixels) and 0.25 
in Vfjoe magnitude as estimated by SEXTRACTOR. This was found 
to produce a total of ~51 000 matched galaxy pairs for each of the 
shear and flexion simulation suites, reducing the galaxy density to 
approximately 45 galaxies arcmin" 2 . More stringent pair exclusion 
criteria (e.g., < 0.09 arcsec separation) were seen to cause rapid re- 
ductions in the numbers of matched pairs. Opting for more tolerant 
criteria produced slightly greater numbers of galaxies (~ 1-3%), but 
did not substantially improve the signal-to-noise of end results. The 
fraction of galaxies lost to catastrophic failures in shapelet mod- 
elling was small (0.4%). 

The difficulty of matching galaxies suggests that signifi- 
cant numbers of objects are being erroneously deblended at the 
SEXTRACTOR detection stage; we found ~ 20-30% losses in this 
simulation, although it should be stressed that deblending errors 
vary and depend on PSF shape, galaxy morphology and SNR. 
Unfortunately, this effect sets a ceiling on the number of objects 
that can be successfully matched after simulating an end-to-end 
pipeline in this way. While including deblending in the simulation 
test contributes realistic noise to the end measurements, it also re- 
duces the constraining power of the simulations as a whole when 
significant numbers of objects are lost, as here. To isolate the ef- 
fects, it may be preferable in future wor k to generate galaxies on 
grids as done in GREAT0 8/GREAT 10 dBridle et alj|2009l l20ld ; 
iKitching et alj|201lll2012h . Flexion and shear estimates were then 
generated for each galaxy in matched pairs, as will now be de- 
scribed. 



6.2 Estimating flexion and shear 

To estimate shear and flexion from shapelet models of galaxies we 
adopt an approach similar to that described by M07, and compare 
the relative values of the shapelet coefficients f n ,m of best-fitting 
shapelet models to derive estimators of lensing distortions. In Ap- 
pendix|A]we describe the generation of flexion and shear estimators 
from catalogues of shapelet coefficient models for a population of 
galaxies with a realistic distribution of fluxes and sizes. We note 
that the multiple-decade variation in these two properties means 
that small modifications to the estimators proposed in M07 must be 
adopted. Described in detail in Appendix lAl we label these estima- 
tors 7, T, and Q for the shear, T flexion and Q flexion, respectively. 

We also note that in order to use these estimators, it was made 
a condition that the shapelet model reach sufficient n max as to con- 
tain non-zero values for all the shapelet coefficients required by 
equations dAlt . dA4t & ( |A9t when estimating 7, T or Q, respec- 
tively. This is to avoid bias from artificially setting coefficient val- 
ues to zero when this is not the most appropriate pri or expectation. 

Overall, this approach differs from that of IVelander et al .1 

d201ll) . who instead take circular profiles and estimate the shear 
and flexion required to distort these objects to m atch the data, in 
a m anner similar to the shear -only estimators of iKuiiker] J2006h 
and iBernstein & Jarvij d2002h . Both methods are similar in that 
they rely to some extent upon shapelet models being an ac- 



curate description of the underlying surface brightness distribu- 
tion of galaxies to avoid what has been i dentified as 'underfil- 
ling bias' in ga laxy shape estimation (e.g. IVoigt& Bridie] l2010l : 
iBernsteinl 12010) . However, in our simulation tests shapelet mod- 
els have been used to provide the underlying galaxy light pro- 
files, and perfect knowledge of the PSF is also available. This 
allows the direct probing of potential biases due to the use of 
noisy data, such as imperfect deblending with SEXTRACTOR, or 
the biased response of non-linear estimators under noise itself (e.g . 
Bernstein & Jarvisll 20021; Iffirata et alj|2004 iRefregier et alj|2012l ; 



Kacprzak et al.ll2012l ; lMelchior & Violall2012h ~ 

Using the method described in Appendix [A] estimators of 
shear and flexion are constructed for each of the galaxies with suc- 
cessful detections and shapelet decompositions described in Sec- 
tion |6. 1 1 In the following Sections, we compare these estimates to 
the input gravitational signal, and explore how the the properties of 
these estimates vary with galaxy properties such as SNR and size. 

6.3 Flexion and shear estimator results 

In Figure[8]we plot (7-7mp.1t) versus y mpM (top panel), (F-Fmpu) 
versus ^input (middle panel), and (Q — Q- mvM ) versus Q mput (bottom 
panel) for our simulated galaxies, with estimators for rotated and 
unrotated pairs of galaxies being me an-averaged to reduce noise 
as described bv lMassev etail d2007bl) . Results were not found to 
differ between components (i.e. results for 71 were consistent with 
those for 72, etc.) and so the plots show both real and imaginary 
parts for each signal. Results were binned according to input signal 
via the image sets described in Section |531 and the points for each 
bin represent the median of the averaged rotated and unrotated es- 
timates in each case. For shear, the median results were consistent 
with results derived from the arithmetic mean; for flexion, the mean 
was found to be very noisy due to the distribution of flexion esti- 
mators (see Section [6~5l l. and so the median was preferred as the 
comparison statistic in both cases. 

The estimated median of each variate in a bivariate distribu- 
tion, taken along each of its two dimensions independently, is a 
good estimator of the central tendency in both variates provided 
that correlations between them are linear. We found no evidence of 
non-linear correlations between the estimates of T\ and J- 2 , indeed 
no evidence of correlations at all, and so do not con sider the use of 
the m edian (rather than, e.g., convex hull stripping; I Velander et al.l 
1201 ll) to be a source of bias in this analysis. The uncertaint ies plot- 
ted sh ow the standard error on the median in each bin dLuptonl 
1 19931) . 

We fit a linear relation to the results of Figure[8] deriving best- 
fitting s lope m and offset c bias parameters as used in the STEP 
project dHevmans et al.120061) : 



7 — 7input = "17input + C, 



(32) 



with similar expressions for the two flexion estimators. The best- 
fitting values, and uncertainties, are given in Figure[8] The real and 
imaginary parts were again found to give consistent results, and 
so the best-fitting m and c describe input versus output for both 
components. We see that the c values are broadly consistent with 
zero for each case, as might be expected given a purely circular 
PSF (although the pixel grid could potentially impart a preferred 
direction also). The offset c was also found to be consistent with 
zero when each i = 1, 2 component was examined independently. 

The value of the multiplicative bias m was found to be 
relatively small in the case of shear, m = 0.028 ± 0.013. 
These results are comparable to those obtained with a number of 
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Figure 8. Lensing measurement calibration results from the full set of simu- 
lated galaxies in matched, rotated pairs, for 7 (top panel), T (central panel), 
and Q (bottom panel). 



shear estimation methods i n the GREAT challenges teridle et al.l 
l201fJ : iKitching et alj|2012n . although there are a number of dif- 
ferences between these simulations and those used in GREAT08 
and GREAT10 (e.g. correlated noise; the distribution of galaxy 
sizes and signal-to-noise; overlapping objects; a purely circular 
PSF). For flexion we detect stronger multiplicative biases, finding 
m — —0.340 ± 0.053 for the median of the T estimators, and 
m = —0.45 ± 0.21 for the median of the Q estimators. Such bi - 
ases are comparable to those identified by I Velander et al. Il l201ll) . 
for galaxies based on analytic profiles and a different distribution 
of sizes and SNR values. We now discuss the variation of our mea- 
sured estimator biases as a function of these properties. 
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Figure 9. Variation of multiplicative bias m (triangular points) and additive 
bias c (diamond points) in shear and flexion estimation versus 'observed' 
SNR for the simulation galaxies (see equation [32}- SNR bins were chosen 
to give equal numbers of galaxy in each bin: the increase in errors for low 
SNR objects is due to the increasing scatter in individual estimates. Solid, 
horizontal lines through points show the extent of each bin. 



6.4 Dependence of bias on noise and apparent galaxy size 

In Figure [9] we plot the dependence of the multiplicative bias m 
and additive bias c (see equation [321 upon observed galaxy signal- 
to-noise ratio (SNR) for pair-matched shear and flexion estimates 
from the simulations. Here, SNR is defined in terms of the SEx- 
TRACTOR output parameters FLUX_AUTO and FLUXERR_AUTO as 



SNR: 



FLUX_AUTO 



X 1/V0.316, 



(33) 



FLUXERR_AUTO 

where the scaling factor 1/V0.316 is taken from lLeauthaud et al.l 
d2007h and adjusts SNR in drizzled HST images to account for cor- 
related noise, adding the assumption that excess Poisson variance 
due to object flux above the background is negligible (it is not in- 
cluded in our simulations). This is only an approxi mate correction, 
based on a simplified model JCasertano et al.l200fj) and an assump- 
tion that the COSMOS drizzling approach closely resembles that 
used in the GEMS noise mosaic (it does, although there is a small 
difference in the kernel used for the Mk II GEMS reduction; see 
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ICaldwell et alj[20o3) . However, it does help take correlated noise 
into account at a level of accuracy that is appropriate given that dif- 
ferences between definitions of SNR can also introduce factor ~2 
discrepancies. 

The SNR ranges for the bins were chosen to give equal num- 
bers of galaxies in each bin; for reference, the bin maxima and 
minima are given in Figure QTJ and indicated in Figure [9] by solid 
horizontal lines. It is difficult to discern strong support for over- 
all trends — results for the lower SNR bins are noisy compared 
to those in higher SNR bins, particularly in the case of the flex- 
ion — although the higher SNR bin consistently gives arguably 
better results. Overal l, flexion results are broadly consistent with 
IVelander et alj d201 ll) . within large errors, hinting that the impact 
of more complex galaxy mor phology (a key differ ence between 
these simulations and those of IVelander et alj|201 lh upon flexion 
estimation bias is not great compared to other properties of the ob- 
servational data. In all cases the additive bias c is consistent with 
zero as would be expected given the circular symmetry of the ap- 
plied PSF 

In addition to the SNR of galaxy objects, the size is an im- 
portant observable quantity for measuring flexion. Unlike shear, 
flexion has dimensions of inverse angle, and larger objects pro- 
vide greater leverage for measuring flexion. In Figure [10] we plot 
m and c for the simulated galaxy sample binned into three bins re- 
flecting angular size, chosen to given equal numbers of galaxies in 
each bin. The size estimate adopted for this binning is V R 2 for the 
best-fitting shapelet models, plotted in arcsec, where the size mea- 
sure R 2 derived from shap elet model coefficients is described by 
iMassev & Refregieil < l2005t) and given in equation iA5l . It is equiv- 
alent to an unweighted, integrated second-moment size over the 
shapelet model. 

Results are qualitatively similar to those found for the SNR 
tests, which is not altogether surprising as we expect a positive cor- 
relation between size and SNR among the simulated galaxy popu- 
lation. There is some evidence that estimators are better for larger 
apparent sizes; the shear in particular now shows a marked trend 
towards improvement for larger objects. But the noise on estimates 
of m for flexion estimators become large for smaller objects, and 
so for flexion clear trends are difficult to discern. 

6.5 The distribution of flexion and shear measurements 

The results of Section |6"4l (Figures |9l & ITOb show a significant in- 
crease in the uncertainty of estimates of m as SNR (and the related 
property V R 2 ) decrease, despite the number of objects being the 
same in each bin. This effect is strongest in the measurements of m 
for T and Q. 

This can be explained by considering the two separable con- 
tributions to measurement uncertainty: the scatter in the intrinsic 
shapes of galaxies prior to lensing, and the noise in shape estimates 
due to noise in pixels. The latter will increase as SNR decreases, 
and may also be a cause o f bias as well as increased uncertainty in 
indiv idual shape estimates jRefregier et al.l2012l ; lMelchior & Violal 
|2012|) . In FigureQTJwe plot the distributions of pair-matched shear 
and flexion estimates as a function of SNR. These paired combina- 
tions cancel the leading order contribution of intrinsic galaxy shape 
to the estimation of each signal, so that the remaining scatter is that 
solely that due to the differing pixel noise applied to each simu- 
lated galaxy. The uncertainty of shear estimates increases as SNR 
decreases, but the effect is far stronger for the estimators of flexion. 

Figure QTJ also illustrates the non-Gaussian distribution of the 
uncertainty in shapelet shear and flexion estimators due to noise 
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Figure 10. Variation of multiplicative bias m (triangular points) and 
additive bias c (diamond points) in shear and flexion estimation ver- 
sus observed size for the simulation galaxies parameterized by V R 2 
iMassev & Re fregiei 20051). I n a similar manner to the SNR binning, size 
bins were chosen to give equal numbers of galaxy in each bin: the increase 
in errors for smaller objects is therefore due to the increasing scatter in in- 
dividual estimates. Solid, horizontal lines through points show the extent of 
each bin. 



in image pixels. For flexion in particular, this distribution is highly 
non-Gaussian, but this is also noticeably true for the shear esti- 
mates: the distributions show a sharp peak in the region of central 
tendency, accompanied by broad wings. However, as the intrinsic 
ellipticity in galaxy images is typically a e ~ 0.3, the shear mea- 
surement uncertainty in the left panel of Figure [TT] only becomes a 
very significant additional contribution for the lowest SNR galax- 
ies. 

For the first flexion T the current best estimate of the un- 
certainty due to intrinsic galaxy shapes is ajr ~ 0.04 arcsec 
dGoldberg & Bacod l2005h . and so in this case the measurement 
uncertainty will dominate over the intrinsic variation over a broad 
range of interest in SNR. The size of the measurement uncertainty 
itself suggests that it may be quite challenging to determine reliable 
es timates of ajr for deeper g alaxy populations than those explored 
in lGold berg & Bacon (2005). The fact that there has been no detec- 
tion of gravitational Q (despite a number of T detections) is also 
plausibly attributable to the extreme measurement scatter seen in 
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Table 4. Best-fitting parameters for the power law model of equation {34), 
describing the NMAD of the shapelet estimator shear and flexion measure- 
ment eiTor as a function of SNR in this study. These models are plotted over 
the data in Figure [12] 
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the right panel of FigureQj}, and while the intrinsic erg is yet to be 
reliably estimated for real galaxies this measurement uncertainty 
is likely to be a significant contribution to noise in future. Flexion 
forecasts which fail to account for the extra uncertainty due to pixel 
noise will provide over-optimistic estimates of future prospects for 
flexion measurement. 

To provide a simple description of the shear and flexion uncer- 
tainty due to pixel noise as a function of SNR, we calculate the Nor- 
malized Median Absolute Deviation (NMAD) for shear and flexion 
measurement error distributions in six bins of SNR. As in earlier 
Sections, the NMAD is useful as a robust estimator of dispersion; 
for distributions such as those in Figure[TT]the population standard 
deviation may not exist, and the sample standard deviation is ex- 
tremely noisy and sensitive to object cuts and outlier removal. We 
plot the NMAD as a function of SNR in Figure [T2] along with the 
best-fitting power law description of these data 

\b 



NMAD = 4io()(SNR/100) 



(34) 



Here SNR =100 has been chosen as the reference scaling value 
since this lies roughly in the middle range of simulated galaxies. 
The best-fitting parameters are give in Table [4] The difference in 
the power law slope between the flexion and shear results, and the 
similarity of the slopes for the T and Q estimators, is interesting: 
it suggests a common origin for the increased flexion measurement 
noise in each case, despite differences in the amplitude of the ef- 
fect. A theoretical explanation for differences between the shear 
and flexion slopes b is unclear (although plausible hypotheses ex- 
ist), but the difference is a clear effect in this data and for these 
estimators. 

These results highlight the importance of considering mea- 
surement noise when discussing flexion estimation in practice. 
We stress again that measurement noise dominates the intrinsic 
a? ~ 0.04 arcsec -1 contribution to uncertainty in T over almost 
the entire range of SNR depicted in Figure[T2](the intrinsic Q is still 
unmeasured). This is a very different regime to that of shear esti- 
mation, where the measurement contribution to uncertainty is dom- 
inated by the intrinsic ellipticity a e ~ 0.3 of galaxies. This is true 
throughout the range of SNR > 10 where shear esti mators have 
been shown have some su ccess in controlling bias jBridle et al.l 
l20irj iKitching et ai]l2012l) . This fact has allowed forecasts for 
shear surveys to proceed using the intrinsic variance alone as a rea- 
sonable approximation to the uncertainty in shear estimators from 
all sources. The results of this study suggest that this should not be 
done for flexion, where noise due to the finite numbers of photons 
arriving at the detector appears to dominate. 



7 DISCUSSION 

We have undertaken a detailed investigation into the problem of 
estimating flexion and shear from noisy galaxy images with sub- 



s'? 1.00 7 



■3 



o 
I 



0.10 : 



0.01 



X I ' 


I 1 


~ s 

N 

s, 


— — y 


X A 


O F - 




-A G 


N> A > 






X : 


-- o "o-. 




o. 


, x : 

"0, 







10 



100 



1000 



SNR 



Figure 12. Normalized median absolute deviation (NMAD) for the distri- 
butions of shear and flexion measurement error in Figure 1111 plotted as a 
function of SNR for the binned galaxy sample. The best-fitting power law 
models are plotted over the data points as straight lines (model parameters 
listed in Tableg). 



stantial variation in underlying galaxy morphology. The simulated 
galaxy models also display a realistic distribution of sizes and ap- 
parent fluxes, drawn directly from the galaxy sample in the HUDF, 
but with additional noise realistic for a wider area (e.g., GEMS; 
COSMOS) survey. 

We have discovered a clear qualitative difference between 
flexion and shear measurement. Whereas noise in shear estimates 
is dominated by intrinsic galaxy ellipticity at typical survey im- 
age depths, the corollary is not true for flexion. Instead, noise at 
the pixel level (due to read noise and finite photon number counts) 
appears to dominate uncertainties in our shapelet flexion estima- 
tors. It will be important to account for this fact when generating 
forecasts of the flexion information content in future surveys: pre- 
dictions based solely on the scatter in intrinsic galaxy flexion ct_f 
alone will be too optimistic. E xisting forecasts such as th ose of B06 
that use aj- ~ 0.04 arcsec _1 dGoldberg & Baconll2005h may only 
correspond to predictions in the limit of large image SNR. 

To aid more realistic forecasts in future we fit a simple power 
law model to the flexion pixel noise as a function of SNR (Sec- 
tion l6.5t . It was found that noise in flexion estimates varies signif- 
icantly more strongly with galaxy SNR than was found for shear 
estimates, nearly a factor of 2 in power law slope. As the majority 
of galaxy images in any wide area survey are likely to be faint, it 
will also be important to consider this effect when forecasting. 

As well as considering the noise in flexion measurements, sys- 
tematic biases in shapelet shear and flexion estimators were also 
investigated. Such tests had only previously been performed with 
galaxy simula tions with less morph ological richness in the under- 
lying models ( tVelander et alj|201lh . The best performing flexion 
estimator in our shapelet pipeline s howed comparable pe rformance 
to the shapelet estimator tested bv lVelander et alj d201ll) . Because 
of the direct use of the HUDF galaxy size and SNR distribution in 
the simulations, exploration of the dependence of biases on these 
properties is less clean, requiring broad bins. We found little strong 
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Figure 11. Distributions of matched measured 71 & 72 (solid and dashed lines, respectively; left panel), T\ & Ti (solid and dashed lines, respectively; centre 
panel), and Q\ & Q 2 (solid and dashed lines, respectively; right panel) from the mean of matched pairs of galaxies in the rotated and unrotated simulations 
(see Section |53V These paired combinations will cancel the leading order contribution of intrinsic galaxy shape to the estimation of each signal: what remains 
will be dominated by uncertainty due to pixel noise. 



evidence of clear systematic trends within the errors for the flex- 
ion results, except that flexion measurement appeared to perform 
somewhat better for large, high-SNR galaxies. Shear results for the 
shapelet pipeline were consistent with the better performing meth- 
ods in the most recent sh ear measurement simulat i on challenges, 
GREAT08 & GREAT 10 jBridle et alj|2009l l20ld ; fetching et ail 
l201lLl2012h . 

In the future, theoretical tools from the increasing body of 
work invested in un derstanding the effects of pixel noise in shear 



measurement (e.g. iBernstein & Jarvis 



Refregier et al .1 12012k 



Kacprzak ct al 



2002); iHirata et alj |2004|; 



20121 ; iMelchior & Violal 



20121 : 1 Miller et alj|2012h will be useful in forming a deeper un- 



derstanding of noise and bias in flexion estimators. In particular, it 
may be possible to understand the apparent qualitative differences 
between flexion and shear estimators under pixel noise. The clear 
difficulties of flexion measurement in practice call for a better un- 
derstanding of such issues if the potential of flexion as a probe of 
small scale power in matter structure is to be realised. 
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APPENDIX A: SHAPELET ESTIMATORS OF FLEXION 
AND SHEAR 

A variety of polar shapelet estimators for shear and flexion are de- 
scribed by M07. In theory, the number of lensing estimators that 
may be constructed using shapelets is only limited by the number of 
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shapelet modes available, n max . However, the majority of shapelet 
models will be highly truncated to somewhat low n max in practice. 
Estimators which make extensive use of higher order information 
prove problematic for many galaxy images, particular those which 
rely upon the convergence (in the sense of converging to a limit) of 
sums over shapelet coefficients (see M07). 

In this paper we limit our investigations of shear and flexion 
measurement primarily to the most simple, and therefore widely 
available in noisy galaxies, shapelet estimators. We now describe 
construction of useful estimators based on the results of M07, tak- 
ing into account important practical realities such as the wide dy- 
namic range of galaxy fluxes and sizes. 



Al Normalizing shapelet coefficients by galaxy size and flux 

The simplest shear estimator is that which employs the polar 
shapelet coefficient / 2i 2 as the galaxy polarization estimator (see 
M07, Section 3, for the precise usage of this term). This quantity 
can in fact be straightforwardly shown to correspond to a Gaussian- 
weighted quadrupole moment of the shapelet model, with a weight- 
ing radius of f3. 

The raw form of the simplest well-motivated estimator that 
can be constructed from /2,2 is as follows: 



~ Gauss /7T 

7 = V2 



(/(),() - /4,o) ' 



(Al) 



where the expectation notation used in the denominator denotes an 
estimat e of the ensemble averag e over all galaxies in the source 
sample dRefregier & Bacorj|2003l : M07). 

However, we have found that this estimator suffers an im- 
portant disadvantage: both the numerator and denominator share 
the dimensionality of the shapelet coefficient s themselves, scaling 
with b oth object flux and inverse scale size /3 dMassev & Refregien 
120051) . This means that estimators of shear (and similarly con- 
structed flexion estimators) vary strongly in magnitude as a func- 
tion of the object size and flux relative to those of the ensemble 
average. 

It is therefore desirable to generalize the estimator of equation 
dAlt to create a new estimator we label 7„, defined as 



7" 



V2 



M/0,0 -/4,o))' 



(A2) 



where the parameter v is suitably chosen to be a property of the 
galaxy image that helps normalize terms in the numerator and de- 
nominator, helping to lessen the impact of the wide dynamic ranges 
of flux and scale length in galaxy samples. This is analogous to the 
adoption of a normalizing factor 1/S in the definition of the im- 
age moments qij, qijk and qijkl in Section [5~3l It should be stated 
from the outset that if v is an estimated property of the image it- 
self it may add both uncertainty and bias to the shear estimator. We 
will discuss how this danger weighs against the benefits of such a 
normalization later in this Section lA3l of this Appendix. 

The simplest possible polar shapelet flexion estimators can be 
constructed from similar combinations of shapelet coefficients. The 
first order T estimator can be expressed as 



^.Gauss 



where 



4 

W (di,i)' 



(A3) 



di,i = \1 - -gjj /o,o + ^2/2,0 + 6/32 £ /a,a - U,o- (A4) 



The shapelet model par ameters e 
iMassev & Refre gier (2005) as 



£(n + l)/„,o, 



and R are defined in 



(A5) 



R 2 F 



where F denotes the shapelet model flux, 

even 

F = /V4^]T/„,o, 



(A6) 



(A7) 



and where all sums are over all even n coefficients / n ,o and / n ,2 
in the model (m = and m = 2 coefficients do not exist for odd 
n). We note that e and R 2 appear only in those terms necessary 
to correctly accoun t for the flexion centroid shift at linear order in 
the applied T (see lGoldberg & Bacon|[2005l : M07). It should also 
be noted that in equation dA4t contains a term proportional to 
e* x fn,2, for n = 2; in M07 these were argued to vanish due to 
rotational symmetry in the source plane. We argue that this term 
should not be omitted since the definition of e above shows the 
inclusion of /2,2 as the first coefficient in the sum. For most galax- 
ies lower-n components of the shapelet model will dominate those 
with higher-n, so that as a first approximation e* x /2.2 ~ |/2,2| 2 - 
This does not cancel due to rotational symmetry across a population 
of galaxies.. 

For the first order Q estimator we have the expression 



~ Gauss 

y 



where 



4VE / 3 , 3 
3/3 (d :i , 3 )' 



(A8) 



d 3 , 3 = (/o,o + /a,o - U,o - /e,o) . (A9) 

Analogously to the case of ^ Gauss ; these estimators of both T 
and Q are directly relate d to Gaussian-weighted octupole moments 
(e.g jQkura et alj|2007l) of the shapelet model galaxy. 

As for shear, the large dynamic range in the shapelet coeffi- 
cients due to varying galaxy size and flux suggests it may be useful 
to rescale numerator and denominator with a generalizing factor v, 
defining: 



T v = 



Q v = 



4 

3/3 (vdi,i) ' 



(A10) 



(All) 



3/3 H, 3 )' 

We now discuss options for a suitable choice of this normalization 
parameter v. 



A2 Choosing a suitable v 

Adopting v = 1 recovers the M07 ^ Gauss estimator. Natural alter- 
native choices for v are combinations of shapelet model parameters 
that make the numerator and denominator of equation JA2t dimen- 
sionless. Two of the simplest and most easily-motivated potential 
combinations are 



v = p/F; 
^ = l//o,o; 



(A12) 
(A13) 
(A14) 



(c.f. the dimen sionless shapelet basis functions introduced in 
iRefregied 12001 .) In Figure I A 1 1 we plot histograms of the values 
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Table Al. Estimates of the skewness T as defined by equation <A15t for the 
distributions of Figure lAll 







Skewness F(x) 




V 


x = v(fo,a - /0.4) 


x = fdi : i 


x = vd it3 


13/F 


0.264 ± 0.002 


0.173 ± 0.002 


0.198 ± 0.001 


V/o,o 


0.0275 ± 0.0002 


-0.0447 ± 0.0004 


0.153 ± 0.001 



of the denominators of the shear and flexion estimators defined by 
equations JA2b . JAlOb & JA lib , from shapelet fits to the simulation 
galaxies in this study. We plot results for the two choices of v ex- 
pressed by equations dA13t & dA12b . The shape of these histograms 
appears to show a consistent trend: using v = l//o,o results in dis- 
tributions of denominator values that have a more clearly-defined 
central tendency than v = ft/F, This fact argues in favour of the 
adoption of l//o,o to normalize to a dimensionless combination of 
shapelet coefficients. 

To explore the distributions further we define the following 
measure of skewness as a combination of the sample mean, arith- 
metic median and standard deviation, commonly known as Pear- 
son's second skewness coefficient: 

r(ar) = 3 [Mean(ar) - Median(x)] /a(x). (A15) 

Estimates of T for the distributions of Figure IA1I with uncer- 
tainties, are given in Table IA1I We see that the skewness in the 
v — f3/F distributions is consistently greater than that in the dis- 
tribution of denominator values when choosing u = l//o,o. Once 
again, this argues in favour of the adoption of v = l//o,o as a 
measurable quantity with which to make the numerators and de- 
nominators of our shear and flexion estimators dimensionless. No 
other combination of observables was seen to provide better per- 
formance, either in terms of showing a clear location for the central 
tendency of the denominator distributions shown in Figure [ATI or 
in terms of reducing the required level of systematic bias calibration 
as shown in Figure [8] To estimate shear and flexion from shapelet 
fits to our simulated images we therefore adopt 7 = 7„, P — T v 
and Q — Q u , where v — l//o,o. 

A3 Potential issues 

The quantity v, which is estimated for each object individually, is a 
function of random variables and therefore a random variable itself. 
Its inclusion in the estimators of equations dA2b . JAlOb & dAl lb 
is therefore a cause of both additional un certainty in estimators o f 
galaxy shape, and potential bias (see, e.g., lMelchior & Violal2012h . 

These undesirable properties must be weighed against the 
practical advantage of the technique in providing dimensionless 
combinations of shapelet coefficients for shape estimation. Given 
the large dynamic range in both apparent galaxy sizes and fluxes in 
extragalactic data, and thus in the raw values of best-fitting shapelet 
coefficients / n , m , this advantage is considerable. Figure lATI shows 
that such combinations can provide distributions with a desirably 
compact support, and the same is true for the denominators of the 
shear and flexion estimators of equations dA2b . dAlOb & dAl It . 

To provide robust estimators of shear and flexion in the pres- 
ence of such a large dynamic range in f n ,m values would probably 
require the splitting of the galaxy sample into cells of apparent siz e 
and flux in a manner similar to that described bv lKaisej | |200C|) . 
Such a scheme carries its own biases, due to the action of shear 
and flexion to carry galaxies between adjacent cells. Given the de- 




Figure Al. Histograms showing distributions of the values averaged in the 
denominators of equations IA21 dAl 01 & dAl It — upper, middle and lower 
panels, respectively — for two choices of normalization parameter v. Solid 
line: choosing v = j3/F. Dot-dashed line: choosing v = l//o,o- 

gree of statistical uncertainty in both shear and flexion measure- 
ments from space-based datasets in the near and medium term, the 
motivation for constructing such dimensional cell-based estimators 
was not strong. However, such a scheme, if correctly implemented, 
might help reduce the calibration factors found in Section [6~3l 
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